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16236 DESIGN CURVES FOR HINGED WAVEMAKERS: THEORY 


KEY WORDS: Curves (geometry); Numerical analysis; Ocean engineering; 
Ocean waves; Theoretical analysis; Theory; Wave measurement; Wave 
propagation; Wave tanks 


ABSTRACT: Dimensionless theoretical design curves for hinged wavemakers of 
variable draft are extended, within the limits of linearized potential wave theory, to 
wave flume geometries consisting of two constant depth domains separated by a 
gradually sloping transition region. These design curves are examined for the following: 
hydrodynamic pressure moment, moment arm, relative moment phase angle, and the 
ratio of the wavemaker stroke to the propagating wave height. The dimensionless 
design curves are shown to demonstrate the effects of wave flume and wavemaker 
geometry on the minimum inertial wave pressure moment on the wavemaker flap 
better than tabular lookup methods. The effects of nondimensional wavemaker 
parameters on the dimensionless wavemaker power curves are explored, and they are 
compared to the energy flux in a propagating periodic wave. Convergence of the 
eigenfunction expansion for the wavemaker velocity for hinged wavemakers of variable 
depth is examined numerically. 


REFERENCE: Hudspeth, Robert T., and Chen, Min-Chu, “Design Curves for Hinged 
Wavemakers: Theory,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. HYS, 
Proc. Paper 16236, May, 1981, pp. 533-552 
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depth; Wave propagation; Wave tanks 


ABSTRACT: Experimental data measured in large-scale wave flume (for the 
dimensionless wavemaker gain function and for the magnitude and phase angle for the 
dimensionless wavemaker hydrodynamic pressure moment on a hinged wavemaker) are 
compared with theoretical values computed from dimensionless design curves which 
are derived from linear wave theory. Dimensionless relative errors between 
experimental and theoretical values are tabulated for 21 discrete wave frequencies in 
each of two water depths. These dimensionless relative errors provide comparisons for 
over one decade of dimensionless relative water depth which cover both intermediate- 
and deep-water wave conditions for wave heights which are approx. 25% of the 
theoretical breaking limit. While measured data from smaller scale wave flumes have 
been previously published for the dimensionless wavemaker gain function, S/H, the 
measured hydrodynamic pressure moment data appear to be unique. 


REFERENCE: Hudspeth, Robert T., Leonard, John W., and Chen, Min-Chu, “Design 
Curves for Hinged Wavemakers: Experiments,” Journal of the Hydraulics Division, 
ASCE, Vol. 107, No. HYS, Proc. Paper 16237, May, 1981, pp. 553-574 


16227 ANALYSIS OF FLOW IN SEDIMENTATION BASINS 


KEY WORDS: Basins (containers); Continuity equation; Finite element 
method; Flow measurement; Flow profiles; Mathematical models; 
Momentum equation; Numerical analysis; Particle distribution; 
Sedimentation; Sedimentation rates; Sedimentation tanks; Turbulent flow; 
Velocity field 


ABSTRACT: The equations governing two-dimensional, steady flow in circular and 
rectangular primary settling basins are formulated. Mean-flow equations expressing 
conservation of mass and momentum are combined with a two-equation turbulence 
model (k-€ model) to compute velocity fields in sedimentation basins. The turbulence 
is characterized by a turbulent kinematic eddy viscosity which is computed directly as 
part of the solution. A numerical solution to the mathematical model is formulated by 
application of the Galerkin finite element method to the equation residuals. The 
resulting non-linear system of equations is solved by a variation of Newton’s method. 
The model is applied to predict the flow pattern in a rectangular settling basin. A 
vector plot of the basin velocity field and contour plot of the turbulent kinematic eddy 
viscosity are presented. 


REFERENCE: Schamber, David R., and Larock, Bruce E., “Numerical Analysis of 
Flow in Sedimentation Basins,” Journal of the Hydraulics Division, ASCE, Vol. 107, 
No. HYS5, Proc. Paper 16227, May, 1981, pp. 575-591 


16245 RESULTS FROM CHANGES IN RIVER FORM 


KEY WORDS: Channel erosion; Drainage basins; Erosion; Flood hydrology; 
Flood plains; Floods; Land use; River flow; Rivers; Sedimentation 


ABSTRACT: Uncertainties resulting from flood-induced major changes in river form in 
some drainage basins are of such magnitude that many problems pertaining to the 
management of water and debris, the management of flood plains, and the effects of 
man’s activities in the basins cannot be resolved except in a probabilistic or general 
way. For affected basins, major changes in channel form can cause significant 
differences in the hydraulics of floods, the surface runoff, recharge to aquifers, erosion 
rates, and sediment yields. Uncertainties develop for several related reasons. The 
capability does not exist to determine adequately whether a river is susceptible to a 
major change in form and to predict when a major river-form change will occur. The 
effects of a local change in river form in a drainage basin usually propagates along 
alluvium-filled valleys to other parts of the basin, but the propagation process is not 
understood. 


REFERENCE: Burkham, Durl E., “Uncertainties Resulting from Changes in River 
Form,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. HYS, Proc. Paper 
16245, May, 1981, pp. 593-610 


16260 CLEAR-WATER SCOUR AND CIRCULAR PIERS 


KEY WORDS: Bridges (piers); Flow profiles; Formulas; Piers (supports); 
Scour; Scouring; Sediment analysis; Streambeds 


ABSTRACT: The potential predictors of the maximum clear-water scour are compared 
with the available experimental data. The range of flow parameters, for which these 
formulas either overpredicted or underpredicted, are delineated. Comparison indicates 
that the scour formula by Laurson and Toch is the best predictor among those 
examined, for it encompasses all data and overpredicts less than the other formulas. 
However, the Laursen and Toch formula predicts that the scour depth is independent 


of sediment size. A new formula for the maximum clear-water scour is proposed; it is 
very similar to that of Laursen and Touch but includes the effect of sediment size on 
scour depth. 


REFERENCE: Jain, Subhash C., “Maximum Clear-Water Scour Around Circular 
Piers,” Journal of the Hydraulics Division, ASCE, Vol. 107, No. HY5, Proc. Paper 
16260, May, 1981, pp. 611-626 
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DESIGN CURVES FOR HINGED 
WAVEMAKERS: THEORY 


By Robert T. Hudspeth," M. ASCE and Min-Chu Chen” 


INTRODUCTION 


The basic wavemaker theory for forced harmonic surface gravity waves derived 
by Havelock (5) has been elaborately extended to design curves by Gilbert, 
et. al (3) for both periodic and random waves for three types of wavemaker 
geometries which are commonly encountered in many modern wave flumes. 
Limited experimental verifications of these design curves were covered by 
Krishnamacker (10) and by Gilbert, et. al (4). Hyun (7) later extended the basic 
theoretical solution for the case of periodic waves only which are generated 
by a hinged-type of wavemaker with a hinge of variable draft. Hyun concluded 
from his theoretical results in the case of deep water waves that the inertial 
component of the hydrodynamic wave pressure force on the wavemaker dominates 
the total pressure force if the hinged wavemaker flap extends over the entire 
fluid depth. Hyun did not present experimental verification of this theoretical 
result; but his conclusion clearly agrees with the design curves previously 
presented by Gilbert, et. al (Ref. 3, Fig. 2, p. 165) for the values of h/gT” 
greater than 0.2. Experimental verifications of these wavemaker solutions for 
various types of wavemaker geometries have been restricted primarily to the 
analytical relationships for the wavemaker gain function which relates wave 
height, H, to the wavemaker stroke, S, which is measured at a particular vertical 
elevation above the bottom [compare to Biesel and Suquet (1) and Ursell, et. 
al (13) for the case of monochromatic linear waves forced by both piston and 
hinged-types of wavemakers; and Madsen (11) and Multer (12) for the nonlinear, 
initial value wavemaker problem formulated by Kennard (8)]. Experimental 
verifications of the theoretical forces, moments, and power requirements for 
wavemakers of any type appear to be relatively scarce. 

The solution presented by Hyun (7) for a hinged-type wavemaker of variable 


‘Assoc. Prof., Dept. of Civ. Engrg. and Ocean Engrg. Programs, Oregon State Univ., 
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FIG. 1.—Definition Sketch for Wave Flume Geometry 


draft is extended in a minor way within the limits of linearized potential wave 
theory to wave flumes which consist of two constant depth regions separated 
by a gradually sloping transition region (Fig. 1). Design curves for the wavemaker 
gain function, S/H, the dimensionless hydrodynamic pressure moment, M’, 
moment arm, //H, and wavemaker power, W’, for this type of wave flume 
geometry are developed which follow the more desirable format presented by 
Gilbert, et. al (3). Experimental verifications of these design curve values were 
measured in a large-scale wave flume at the Oregon State University Wave 
Research Facility (OSU-WRF) and are presented in the companion paper (6). 


Wavemaxer THEORY 


For convenience in developing universal design curves, all variables will be 
made dimensionless at the outset by the usual physical variables; viz., a length 
scale = g/w’; a time scale = w~', and a mass density scale = p, in which 
g = gravitational constant; w = radian frequency; and p = mass density of 
fluid in the wave flume. The relationships between physical variables (denoted 
by *) and dimensionless variables are given by the following: 

2 


(x,y,S,h,H) = (x*,y*,S*,h*, H*) — 


(K,,,«) = (Ki ,«*)—> 


(t,T) =(t*,T*)o 


7 
w 


F =(F*)—~; M=(M*)—~; (W) = ((W*)) — (le) 
Pg pg 


Pg 

The two-dimensional, irrotational motion of an incompressible, inviscid fluid 
in the wave flume geometry shown in Fig. | for a right-handed Cartesian coordinate 


system may be obtained from the directional derivatives of a dimensionless 
scalar velocity potential, ®(x, y,¢), according to 


u(x,y,t)= -®, 
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n 
i /, An 
| 
g 
(15) 
| (1c) 
© = (6*) — ) 
| g 
(20) 
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v(x,y,t) = 
in which the subscripts denote partial differentiation with respect to the indepen- 


dent spatial variables denoted by the subscript. The total pressure may be 
determined from the linearized Bernoulli equation given by 


P(x, y,th=®,-y 


For simple harmonic wavemaker motions which are strictly periodic in period 
T = 2, a dimensionless spatial velocity potential, (x,y), may be defined 
by the real part of 


P(x, y,t) = d(x, y) exp (it) 


which must be a solution to the dimensionless equation of continuity given 
by 


o,.,+6,,=0; O=x<+0, -h=y=0 
Considering first the constant depth region immediately adjacent to the 


wavemaker in Fig. 1, the linearized homogeneous boundary conditions for the 
constant horizontal boundaries are given by 


y=0,x>0 

6,=0; y=-h,x>0 

Eqs. 5 and 6 form a well-posed Sturm-Liouville problem for @ by the method 
of separation of variables. The multiplicative coefficients in the solution for 
this type of eigenvalue problem must be quantified by an inhomogeneous bound- 


ary condition which is given by the prescribed motion of the wavemaker. 
Expressing the wavemaker surface by the following Stokes’ material coordinate: 


S=x-x(y,t)=0 


the inhomogeneous boundary condition for the fluid motion may be obtained 
from the Stokes’ material derivative of Eq. 7 according to 


Expanding Eq. 8 in a Maclaurin series about x = 0 and retaining only the 
linear terms yields the following linearized inhomogeneous boundary condition 
for small oscillations of the wavemaker: 


= x =0; 


A radiation condition that requires only outgoing, right-progressing waves in 
Fig. 1 completes the linearized boundary value problem. If the reflections from 
both the gradually sloping transition region between the two constant depth 
domains and from the sloping terminal beach are neglected, the conservation 
of energy flux [compare to Eagleson and Dean (2)] may be used to compute 
the fluid motion in the second constant depth domain of depth, A(1 — A). 

Specifically, the simple harmonic displacement of the hinged wavemaker in 
Fig. 1 may be given by the real part of 


x(y, = &(y) exp i(t + a) 


=0 
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in which a = an arbitrary initial phase angle; and 


Ss y 
&(y) <( ( (lla) 


The arbitrary initial phase angle, a, has been introduced in anticipation of data 
analyses by a finite Fourier transform (FFT) algorithm in order to more efficiently 
determine the phase and magnitude of the measured simple harmonic physical 
variables. Restricting the phase angle, a, to integer multiples of nm/2 radians 
as in the problem formulations in other references is too restrictive and 
unnecessarily complicates data analyses by FFT methods. 


The well-known (1) velocity potential solution to Eq. 5 which satisfies Eq. 
6 is given by the real part of 


n=l 
in which the unknown multiplicative coefficients A, are either pure real or 


pure imaginary. The orthonormal eigenfunctions /,(y) are given by the following 
set [Kreisel (9)]: 


cos K,(y +h) 


in which the normalizing constants, N,,, are determined from 


2K,h+sin2K,h 

N,= cos’ K,(y + hk) 1,2, .... AS) 
4K, 

provided that the eigenvalues K,,h are computed from 


The orthogonality of the orthonormal eigenfunctions f, in the interval of 
orthogonality —h =< y =< 0 allows the multiplicative unknown coefficients A, 
to be computed in a best least-square sense from the following integral: 

0 


which yields for the geometry given in Fig. 1 


S  [K,h(1 -—8)sin K,h + cos K,h — cos K, hd] (:=) 
(K,)°(2K,,h + sin2K,h)'? 
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Again, the phase shift 7/2 has been introduced in order to reference the velocity 
potential to the displacement of the wavemaker in anticipation of data analyses 
by an FFT algorithm which is referenced to the displacement of the wavemaker 
instead of to the velocity of the wavemaker as prescribed by the boundary 
value problem in Eq. 9. 

The simple harmonic linear free surface profile, (x,t), of amplitude H/2, 
may be determined from the linearized Bernoulli equation given by Eq. 3 within 
the first constant depth region, h, immediately adjacent to, but sufficiently 
far away from, the wavemaker in order for the evanescent modes (i.e., nm = 
2) to have decayed to less than 1% of their values at the wavemaker, (x > 
3h, say); i.e.: 


which gives A, =——— 


Substitution of Eq. 18 into Eq. 20 may be shown to reduce to the gain function 
ratio H/S, obtained by Hyun (Ref. 7, Eq. 15, p. 3) for the special case that 
the stroke of the hinged wavemaker, S, is measured at { = 0, as well as to 
the ratio obtained by Ursell, et. al (Ref. 13, Eq. 3.10, p. 37) for the special 
case that 8 = ¢ = 0. However, many wave flume experiments are designed 
to be conducted in the second domain of constant depth A(1 — A) in Fig. 


1. Invoking the conservation of energy flux from linear wave theory (2) yields 
for the wavemaker gain function 


2«h(1 — A) 
+ 
KhH |* a sinh 2xh(1 — A) 
2Kh 


4sinh Kh Lx 


sinh 2 Kh 
2 Kh + sinh 2Kh 
Kh(1 — 8) sinh Kh — cosh Kh + cosh Khid 


in which the propagating wave number, «, in the constant depth region A(1 
— A) is determined from 


— A) = — A) tanh xh(1 — A) 


Fig. 2 demonstrates this design wavemaker gain function for various wave 
flume geometries; and demonstrates that the effect of reducing the wavemaker 
draft (i.e., increasing 5) results in increasing the magnitude of the wavemaker 
gain function. This increase is less for the deep water case and supports a 
similar conclusion reached by Hyun (7). 

The moment about the wavemaker hinge which is induced on the wavemaker 
of width B shown in Fig. 1 by the fluid side of the wave flume only is simply 
obtained from the hydrodynamic pressure moment computed from the linearized 
Bernoulli equation. For a positive moment vector defined by M = r x F, 
the component of torque induced on the wavemaker by the fluid will yield 
a positive moment defined by the real part of 


| 
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FIG. 2.—Dimensioniess Wavemaker Gain Function, S/H, as Function of Wave Flume 
Geometry 


| [y + A(1 —8)] y) dy 


1) 


M = —M, sin (t + a) — M, cos (t + a) 
M = |M| cos (t + ayy) 


in which a,, = measured phase angle for the moment. The amplitude of the 
moment due to the propagating mode is given by 


2x h(1 — A) 
1/2 1 + 
sinh 2x h(1 — A) 
sinh Kh 2Kh 
sinh 2Kh 
[Kh(1 — 8) sinh Kh — cosh Kh + cosh Kh8] 
and the amplitude of the moment due to the evanescent modes is given by 
— A) 
1/2 + 
HBK I<] sinh 2xh(1 — A) 
2sinh Kh 2Kh 
sinh 2Kh 
(2Kh + sinh 2Kh) 
[Kh(1 — 8) sinh Kh — cosh Kh + cosh Kh8] 


K 


K 


538 HY5 
10.00 
"100 
| [1.0 11.0 [1.0 [1.0 [1.0 [1.0 110 [1.0/1.5 [45 [15 05 [as los 2 
8 | a0 |00 | .25 | 25 | 25 0.5/0.5 |05 [0.0/0.0 |0.0 |0.0 [0.0 | L 
° 
10 1.00 
Q= 
Lo 


HY5 HINGED WAVEMAKERS: THEORY 


> [K,A(1 —8) sin K,h + cos K,h — cos K, hd]? 
K4(2K,h + sin2K,h) 


Eqs. 24a and 24b may be shown to reduce to Eqs. 20b and 20c, respectively, 
given by Hyun (Ref. 7, p. 4) for the special case that the wave flume has 
a constant depth, that the wavemaker stroke is measured at the still water 
level, and finally, that fluid exists on only one side of the wavemaker. 

The relationships between the magnitude of the measured hydrodynamic 
pressure moment, |M|, and the magnitudes of the theoretical components of 
the hydrodynamic pressure moment on the wavemaker are easily shown to 
be the following (compare Appendix I): 


M, = |M| sin (a, — @) 
M,=—|M| cos (ay, — @) 


The theoretical relative moment phase angle, B,,, between the hydrodynamic 
moment and the displacement of the wavemaker is given by 


B,, = arctan M. 


which may be computed from the measured FFT coefficients which are referenced 
to the arbitrary initial displacement of the wavemaker by subtracting the two 
measured phase angles from each measured signal according to (compare to 
Appendix I) 


Bu = Oy 


The components of the rotational velocity of the wavemaker about the hinged 
bottom may be determined from 


(r x q)-e, = —y(u, + u,) 


By substituting for the components of the horizontal velocity from Eqs. 2a, 
13, and 18 and by equating real parts, it may be shown that the propagating 
component of the wavemaker moment, —M, sin (f + a) given by Eq. 235 
is in-phase with the wavemaker angular velocity; while the evanescent component 
of the wavemaker moment, —M, cos (t + a), is in-phase with the wavemaker 
angular acceleration. 

Similarly, the total hydrodynamic pressure force in the x direction which 
is induced on the wavemaker by the fluid side only of the wave flume may 
be computed by again integrating over the wetted area of the wavemaker the 
dynamic component of the linearized Bernoulli equation; i.e., from the real 
part of 

0 
o(0,y) dy 
A(é—1) 
F = F, sin (¢ + a) + F, cos (t + a) 


F = |F| cos (t + a,) 


in which the a, = measured phase angle for the pressure force. The amplitude 
of the dynamic pressure force due to the propagating mode is given by 
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4) 
1/2 1+ 
HB K sinh 2xh(1 —A) 
|— (sinh Kh sinh ) (27a) 
2K’ sinh Kh Lx 2Kh 
+ 
sinh 2Kh 


and the amplitude of the dynamic pressure force due to the evanescent modes 
is given by 


2xh(1 — A) 
+ 
HBK sinh 2x A(1 — A) 
F, = ——— | — 
2 sinh Kh L x 2Kh 
sinh 2 Kh 


(2Kh + sinh 2Kh) 
[Kh(1 —8) sinh Kh — cosh Kh + cosh Khd } 


[K,,A(1 — 8) sin K,h + cos K,h — cos K,, hd] 


(sin K,h — sin K,hd) (27b) 
K,(2K,h + sin2K,h) 


Eqs. 27a and 27b may again be shown to reduce to Eqs. 216 and 2lc, 
respectively, given by Hyun (Ref. 7, p. 4) for the previously noted special 
geometry. In addition, Eq. 27a may be shown to reduce to the dimensionless 
propagating force component (termed resistive component) for the special constant 
depth case treated by Gilbert, et. al (Ref. 3, Table 1, p. 184), but not to their 
dimensionless evanescent pressure force component (termed inertial component). 
The reason for this lack of equivalence for the evanescent modes is not known. 
Numerical values estimated from their design curve (Ref. 3, Fig. 2, p. 165) 
do, however, agree with numerical values computed from Eq. 276 for the 
equivalent special wave flume geometry. 

The relationships between the magnitude of the measured hydrodynamic 
pressure force, |F|, and the magnitudes of the theoretical components of the 


total hydrodynamic pressure force are easily shown to be the following (compare 
to Appendix I): 


F, = |F| sin (a, — a) 


The theoretical relative force phase angle, B,, between the total hydrodynamic 
pressure force and the displacement of the wavemaker is given by 


F, 


which may be computed from the measured FFT coefficients which are referenced 


to the displacement of the wavemaker by subtracting the two measured phase 
angles according to (compare to appendix I) 
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wavemaker velocity computed from u = u, + u,, it may be shown that the 
propagating component of the wavemaker force, F, sin (¢ + a), given by Eq. 
26b is in-phase with the horizontal component of the wavemaker velocity; while 
the evanescent component of the wavemaker force, F, cos (t + a), is in-phase 
with the horizontal component of the wavemaker acceleration. Note that these 
are the horizontal components of the angular velocity of the wavemaker. 

The amplitude values computed from the evanescent eigenmode summations 
in Eq. 27b become negative for summations in relative water depths less than 
h(i — A)/L, < 0.4. This indicates that for most intermediate and shallow 
water conditions, the hydrodynamic pressure force leads the wavemaker velocity 
in a phase vector diagram; while for deep water relative depth conditions, the 
hydrodynamic pressure force lags the wavemaker velocity phasor. These phasor 
relationships are shown schematically in Fig. 3. Summation values for the 
evanescent eigenmode component of the moment never become negative over 


Fw) 
[2 > 04] 


FIG. 3.—Phase Vector Diagram for Hinged Wavemaker Pressure Force for Wavemaker 
Displacement Phase Vector S(t) = S/2 exp i(wt + a) 


the two decades of dimensionless frequency which were computed. We note 
that the relative phase angles computed by Eqs. 25c and 28c are referenced 
to the wavemaker displacement and are, therefore, complementary phase angles 
to those given by Hyun (Ref. 7, Table 2, p. 6) which are referenced to the 
wavemaker velocity. 

Finally, neglecting the viscous and mechanical losses of the wavemaker system, 
the average rate of power required by the wavemaker to do work on a column 

f fluid on only one side of the wavemaker may be computed by temporally 
averaging the rate of work done by the hydrodynamic pressure force on a 
column of fluid according to 


(W)=B | (R [® ,(0,y, t)] R[-® ,(0,y,t)]) dy 
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in which the brackets, (-), denote the temporal averaging operator for simple 
harmonic motion and R denotes the real part of a complex-valued quantity. 
Due to the orthogonality of the orthonormal eigenfunctions, f,(y), and to the 
periodicity of the averaging operator, it can be shown that the evanescent 
eigenmodes do not contribute to the average rate of work and that 


BH’ 2xh(1 — A) 
16 Kx tanh Kh sinh 2xh(1 — A) 
Desicn Curves ror Hinceo Wavemakers 


The dimensional design parameters which are most frequently specified for 
experimental work in wave flumes are the design wave height, H, design water 
depth, Ah, and design wave period, 7. Consequently, it is more desirable to 
nondimensionalize the dimensionless wavemaker variables which were previously 
derived by these three design parameters instead of by the previously given 
universal parameters which are independent of the design parameters of wave 
height and wavemaker stroke. Dimensionless design curves for the wavemaker 
variables which have been nondimensionalized by these three design parameters 
may be used more directly in wave flume experimental work. Note that we 


TABLE 1.—Dimensionless Parameters for Hinged Wavemaker (Repeating Variables: 
Fluid Density, p; Gravitational Constant, g; and Wavemaker Fluid Depth, A) 


(7) 


W/ogh? (gh)'” H/h | 


do not normalize the wavemaker variables by the wavemaker stroke, S, since 
this parameter is not known a priori but must be computed from the wavemaker 
gain function using the specified design wave height. The numerical computations 
which are needed to construct the dimensionless design curves require summing 
an infinite series for the evanescent eigenmodes. These series summations were 
truncated when the addition of a single term failed to change the sum by more 
than 0.1% (compare to Gilbert, et. al (3)]. We also return to dimensional variables 
but omit the superscript asterisk (*) for convenience. 

Table | summarizes the dimensionless parameters obtained from a dimensional 
analysis for forced harmonic surface gravity waves generated in the wave flume 
shown in Fig. | if the repeating variables are chosen to be the fluid density, 
p, the gravitational constant, g, and the wavemaker fluid depth, h. Forming 
various ratios with these dimensionless parameters yields the following dimen- 
sionless wavemaker design curve variables: 


M 


l 
—8) 


(1) (3) 


HINGED WAVEMAKERS: THEORY 
F 


1 
(W) 


1 
Bigh)'” 


in which the superscript primes denote dimensionless variables and the deep 
water wave length computed from linear wave theory is 


The amplitude of the dimensionless hydrodynamic moment on the hinged 
wavemaker due to the propagating mode [termed resistive component by Gilbert, 
et. al (3)] is now given by 


[Kh(1 — 8) sinh Kh — cosh Kh + cosh Khd] [E !’ sinh 2xh(1 — A) 


= 200 (33a) 


(Kh)°(1 — 8)(1 — A) sinh Kh 2Kh 
+ 
sinh 2Kh 


K 


and due to the evanescent modes [termed inertia component by Gilbert, et. 
al (3)] by 


M! = 200 


(Kh)(2Kh + sinh 2Kh) 
(1 — 8)(1 — A) sinh Kh [Kh(1 —8) sinh Kh — cosh Kh + cosh Kh8] 


K 


sinh 2kA(1—A)| [K,A(1—8)sin K,h + cos K,h — cos 
2Kh (K,h)*(2K,h + sin 2K,h) 
, sinh 2Kh 


(33d) 


The magnitudes of the dimensionless evanescent, propagating and total hydro- 
dynamic moments are shown in Fig. 4 for various wave flume geometries. 
By nondimensionalizing the hydrodynamic moment on the wavemaker by the 
three previously identified design wave parameters of H, 7, and h, Eqs. 33a 
and 335 now become independent of the dimensionless height of the wavemaker 
stroke above the hinge, H. Note also that the minimum value of the dimensionless 
hydrodynamic moment is easily located graphically in this Gilbert, et. al (3) 
type of figure in contrast to the diagonal term look-up method presented by 
Hyun (Ref. 7, Table 2). 

The amplitude of the dimensionless hydrodynamic pressure force due to the 
propagating mode [termed resistive by Gilbert, el. at (3)] is now given by 


h(i — A) 
(31d) 
L, (32) 
> 
| 
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2«xh(1 — A) 
(sinh Kh — sinh Khd ) E , sinh 2xh(1 —A) 
(Kh)?(1 —8)(1 — A) sinh Kh 


(34a) 


2Kh 
+ 
sinh 2Kh 


and due to the evanescent eigenmodes [termed inertia by Gilbert, et. al (3)] 
by 


2nM Kh(2Kh + sinh 2Kh) 
* = 8)(1 A) sinh KA[Kh(1 —8) sinh Kh — cosh Kh + cosh Kh8] | 


2xh(1—A) 
sinh 2xh(1 — A) s [K,h(1 —8)sin K,h + cos K,h — cos K, hd] 


2Kh (K,h)°(2K,h + sin 2K,h) 
n 


1 + 
sinh 2 Kh 


(sin K,A — sin K,, hd) 


Hyun (Ref. 7, Table 2, p. 6) tabulated the relative contributions from each 
of the two components of the total hydrodynamic pressure force and moment 
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FIG. 4.—Dimensioniess Hydrodynamic Moment, M’, as Function of Wave Flume 
Geumetry 


as a function of the dimensionless relative water depth. These relative effects 
were determined from the phase angles between the hydrodynamic pressure 
force or moment and the velocity of the hinged wavemaker. Since the two 
components of pressure due to the propagating (or resistive) and evanescent 
(or inertia) modes are 90° out-of-phase with respect to each other, these 
comparisons of relative phase angles reflect the relative importance of each 
hydrodynamic component to the total pressure force or moment. This comparison 
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FIG. 5.—Relative Phase Angle, 8,,, between Total Hydrodynamic Moment and 
Displacement of Wavemaker 
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FIG. 6.—Relative Phase Angle, B ., between Total Hydrodynamic Force and Displace- 
ment of Wavemaker 
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is more easily demonstrated in the Gilbert, et. al (3) type of graph in which 
the relative theoretical phase angle may be computed from Eq. 25c for the 
theoretical hydrodynamic moment for the special case of zero arbitrary initial 
phase angle or a = 0. Figs. 5 and 6 demonstrate the relative phase angle for 
the hydrodynamic pressure moment and force, respectively in a Gilbert, et. 
al (3) type of graph for various wave flume geometries. The relative importance 
of the evanescent eigenmodes for deep water conditions as a function of 
wavemaker geometry is readily apparent by comparing these phase angles with 
the phasor diagram in Fig. 3. 

In the design curves presented by Gilbert, et. al (3), the moment arm for 
the hydrodynamic pressure force, /, was assumed to be equal to the wave 
flume water depth at the wavemaker. A theoretical dimensionless moment arm 
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FIG. 7.—Dimensionless Hydrodynamic Moment Arm, //H, as Function of Wave Flume 
Geometry 


may be estimated from a ratio of the maximum amplitude of the hydrodynamic 


moment to the magnitude of the hydrodynamic force at the same instant of 
time according to 


in which the magnitude of the hydrodynamic force is evaluated at t = a,, 
(see Eq. 26c). 

Fig. 7 demonstrates the effects of the wavemaker and the wave flume geometries 
on the dimensionless theoretical moment arm. The ratio of hydrodynamic moment 
arm to wave flume water depth is less than unity for the wave flume geometry 
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analyzed by Gilbert, et. al (3) (viz, H = 1.0, 8 = A = 0.0) which indicates 
that, theoretically, their values are conservative for design applications. 

Hyun (Ref. 7, Fig. 3, p. 5) also presents a graph of the dimensionless average 
wavemaker power versus a dimensionless relative water depth and concludes 
that this dimensionless quantity is independent of the wavemaker geometry 
(viz, the height of the wavemaker hinge and the location of the measurement 
of the wavemaker stroke above the bottom). It is well-known that if the frictional 
and mechanical losses in a wavemaker are neglected, then the average rate 
of power required by a wavemaker is identically equal to the average rate 
of energy flux in the propagating wave component. If the average rate of power 
required by any wavemaker is nondimensionalized by the average rate of energy 
flux in the propagating wave, then this dimensionless quantity is a constant 
equal to unity which is obviously independent of the wavemaker geometry. 


1.00 


FIG. 8.—Dimensionless Average Wavemaker Power, (W’ ), as Function of Wave Flume 
Geometry 


From a dimensional analysis, it may be readily shown that a dimensionless 
average wavemaker power may be computed from 


2xh(1 — A) 
(W’) =| —— | 1 + 
2(1- A) sinh 2xh(1 — A) 


Fig. 8 demonstrates the effects of the wave flume geometry on the dimensionless 
wavemaker power. The special case of A = 0 corresponds to the wave flume 
geometry treated by Hyun (7) provided his dimensionless power is divided by 
the square root of the gravitational constant, g. Evidently, the choice of variables 
used to nondimensionalize the dimensionless average wavemaker power is 
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responsible for the functional dependency of this design variable on the geometry 
of the wavemaker. 


Conc.usions 


The design curves presented by Gilbert, et. al (3) for the generation of periodic 
waves by hinged wavemakers are extended to variable-draft hinged wavemakers 
in wave flumes consisting of two constant depth domains which are connected 
by a gradually sloping transition region in which wave reflections are assumed 
to be negligible. These types of design curves when nondimensionalized by 
the design parameters of H, T, and h, are found to present the important design 
considerations better than the tabular method presented by Hyun (7). 

Specifically, the dimensionless wavemaker gain function, S/H, shown in Fig. 
2 for various geometries is shown to increase with an increase in the dimensionless 
hinge height, 5, for a fixed value of wavemaker stroke height, H. The effect 
is simply a vertical displacement upward of this dimensionless design curve 
without change in shape. On the other hand, if the dimensionless hinge height, 
8, is increased while the wavemaker stroke is held at a constant elevation (i.e., 
H decreased), the dimensionless wavemaker gain function design curve is 
decreased without change in shape. 

The dimensionless relative water depth which locates the minimum value 
for the dimensionless moment amplitude due to the evanescent eigenmodes, 
M! [termed interia by Gilbert, et. al (3)], is shown in Fig. 4 to increase with 
increasing dimensionless hinge height, 5, and increasing constant dimensionless 
transition height, A. The dimensionless maximum total moment amplitude, M’, 
is almost totally due to the propagating eigenmode, M? , for intermediate water 
depth and almost totally due to the evanescent eigenmodes, M/ , for deep water 
depths. The minimum dimensionless total moment is easily found from this 
Gilbert, et. al (3) type of design curve. 

The design curves for the theoretical relative phase angles for the hydrodynamic 
pressure moment, B,,, and for the hydrodynamic pressure force, B,, shown 
in Figs. 5 and 6, respectively, support the relative importance demonstrated 
by Hyun (7) of the amplitude of the propagating mode for intermediate water 
depths (i.e., 2 < 0.5 and B ~ 90°) while the evanescent eigenmodes dominate 
for deep water conditions {i.e., 2 > 0.5 and B ~ 0°}. These phase angles 
are referenced to the displacement of the wavemaker and are, therefore, 
complementary phase-angles to those computed by Hyun (7). 

The parametric dependency of the dimensionless power, W’, for a hinged 
wavemaker on the geometry of the wave flume is shown to depend upon the 
variables used to nondimensionalize the wavemaker power and caution should 
be exercised when drawing conclusions about the power requirements for hinged 
wavemakers of variable draft [compare Hyun (7) p. 4 and Fig. 3, p. 5]. 

Experimental verifications of these theoretical dimensionless design curves 
were obtained in the OSU-WRF for wave heights which are approx 25% of 
the theoretical breaking wave height and the results are presented in the companion 
paper (6). 

Since it is the dimensionless wavemaker moment, M’, which is usually required 
for the design of hinged wavemakers and since experimental verifications in 
the OSU-WRF were obtained for the wavemaker moment in order to avoid 


| 
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having to estimate the moment arm for the hydrodynamic pressure force, no 
design curves for wavemaker forces have been presented. The amplitudes for 
the evanescent eigenmode forces given by Eq. 27b become negative for relative 
water depths less than h(1 — A)/L, < 0.4; while no negative summations were 
obtained for the evanescent component of the moment. In order to obtain an 
accuracy of 0.1%, 27 evanescent eigenmode terms were required for the most 
extreme deep water case computed by Hyun (Ref. 7, Table 2, p. 6). Most 
other cases in intermediate water depths required less than five evanescent 
terms in order to obtain 0.1% accuracy. 
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Aprennix 1.—Reso.ution oF THeoretica, PHase Ancies From FFT Data 


The theoretical relative phase angle between the displacement of the wavemaker 
and the hydrodynamic pressure moment may be readily computed by first 
expanding the hydrodynamic pressure moment given by Eq. 23c according to 
the following angle sum identity: 


|M| cos (t + ay) = |M|costcosa, — |M|sint sina, 


Expanding Eq. 236 similarly and equating coefficients of sin t and cos t separately 
yields the following two equations: 


—M, cosa + M, sina = —|M| sina 
—M, sina — M, cosa = |M| cosa y 


which gives the following solutions for the amplitudes of the components of 
the hydrodynamic moment: 


M, = |M| sin (a,, — a) = |M| sin (B,,) 
—M,=|M| cos (ay, — «) = |M| cos (By) 


in which the theoretical phase angle, B ,,, is relative to the wavemaker displacement 
and Eq. 25d follows. Resolution of the theoretical relative phase angle for the 
hydrodynamic pressure force on the wavemaker, 8, follows in an analogous 
manner. 
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Aprennix 


The following symbols are used in this paper: 


A, multiplicative coefficient for mth eigenfunction of 
velocity potential; 
total width of wavemaker; 
unit vector in right-handed Cartesian coordinate sys- 
tem; 
dimensionless total hydrodynamic pressure force; 
hydrodynamic pressure force amplitude due to 
evanescent eigenmodes; 
hydrodynamic pressure force amplitude due to propa- 
gating eigenmode; 
orthonormal eigenfunction; 
gravitational constant; 
deterministic wave height; 
still water depth at wavemaker; 
imaginary unit; 
wave number for propagating eigenmode in constant 
depth domain adjacent to wavemaker (= iK , ); 
wave number for evanescent eigenmodes (n = 2); 
deep water wave length; 
moment arm for total hydrodynamic pressure force 
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measured vertically above wave flume bottom; 
dimensionless total hydrodynamic moment on wave- 
maker; 
dimensionless hydrodynamic moment amplitude due 
to evanescent eigenmodes; 
dimensionless hydrodynamic moment amplitude due 
to propagating eigenmode; 
normalizing constant for mth orthonormal eigenfunc- 
tion; 
summation index for orthonormal eigenfunctions; 

Pp hydrodynamic pressure; 

q = ue, + ve, + we, velocity vector; 
r= xe, + ye, + ze, position vector; 

wavemaker stroke measured at height of wavemaker 
piston; 
Stokes’ material coordinate for wavemaker displace- 
ment; 
wave period; 
time; 
horizontal component of water particle velocity; 
vertical component of water particle velocity; 
average wavemaker power; 
horizontal component of water particle velocity nor- 
mal to x-y plane; 
horizontal coordinate axis with origin located at 
undisturbed water level at wavemaker; 
vertical coordinate axis with origin located at still 
water level at wavemaker; 
horizontal coordinate axis which is normal to x-y plane 
in right-handed Cartesian coordinate system; 
arbitrary initial phase angle for wavemaker displace- 
ment due to FFT digitizing process; 
measured relative phase angle for measured hydro- 
dynamic pressure force; 
measured relative phase angle for measured hydro- 
dynamic pressure moment; 
theoretical relative phase angle between hydrody- 
namic force and displacement of wavemaker; 
theoretical relative phase angle between hydrody- 
namic pressure moment and displacement of wave- 
maker; 
dimensionless vertical transition height between two 
constant depth wave flume domains; 

5) dimensionless height of wavemaker hinge above bot- 
tom; 

4 dimensionless depth below still water level to wave- 
maker piston; 

H =(1-8- 9 dimensionless height of wavemaker piston measured 

above wavemaker hinge; 
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instantaneous water surface elevation measured posi- 
tive upwards from undisturbed still water level; 
wave number for propagating eigenmode in constant 
depth test domain; 
vertical dependence of prescribed wavemaker dis- 
placement; 
fluid density; 
temporal and spatial scalar velocity potential; 
spatial scalar velocity potential; 
prescribed wavemaker displacement; 
relative water depth; 
radian frequency; 
V = (a/ax)e, 
+ (d/dy)e, 
+(d/dz)e, = gradient operator; and 
( ) averaging operator. 


Superscripts 
dimensional variable; and 
dimensionless dependent variables. 


Subscripts 

evanescent mode (n = 2); 
hydrodynamic pressure force; 
hydrodynamic pressure moment; 
eigenfunction mode summation; 
deep water conditions; 
propagating eigenmode (n = 1); 
x coordinate axis; 

y coordinate axis; and 

z coordinate axis. 
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DESIGN CURVES FOR HINGED WAVEMAKERS: 
EXPERIMENTS 


By Robert T. Hudspeth," M. ASCE, John W. Leonard,’ A. M. ASCE, 
and Min-Chu Chen* 


INTRODUCTION 


The theory for the elaborate design curves for hinged wavemakers presented 
by Gilbert, et. al (4) has been extended within the limits of linearized potential 
wave theory for hinged wavemakers of variable draft in the companion paper 
(6) for wave flume geometries which consist of two horizontal constant depth 
domains connected by a gradually sloping transition region (compare to Fig. 
1). These Gilbert, et. al type of design curves were shown in the companion 
paper to demonstrate the effects of wavemaker and wave flume geometry on 
the design curve variables better than the tabular look-up method presented 
by Hyun (8) for hinged wavemakers of variable draft. However, experimental 
verifications of these design curve variables for variables other than the dimen- 
sionless wavemaker gain function, S/H, appear to be scarce. 

The design curve variables for the wavemaker gain function, S/H, dimension- 
less hydrodynamic moment magnitude, M’, and relative phase angle, B,,, were 
measured in the Oregon State University—-Wave Research Facility (OSU-WRF) 
in two water depths for over one decade of dimensionless wave frequency. 
Comparisons with theoretical values for those wavemaker design variables which 
are presented in the companion paper indicate generally good agreement for 
deterministic waves which are approx 25% of the theoretical breaking wave 
height. These comparisons of measured data with the magnitudes and phases 
of the hydrodynamic moment on a hinged wavemaker are believed to be unique. 


Wave Fiume Description 


The theoretical design curve values derived in the companion paper are 
compared with a unique set of data recorded in the OSU-WRF shown in Fig. 
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1. The wave flume is 89.64 m tong, 3.66 m wide, and has a relocatable bottom 
in the test section. Experimental data were obtained for two water depths at 
the wavemaker of h = 3.96 m and 4.42 m. 

The hinged wavemaker in the OSU-WRF is sealed under pressure along both 
vertical sides and the horizontal hinged bottom so that water is not required 
in the dry well behind the wavemaker. This dry well reduces the power 
requirements for the wavemaker and also eliminates the need to place any 
wave dissipating material behind the wavemaker. The hinged wavemaker is 
controlled by a 112-kW, 20.7-MPa hydraulic pump through a hydraulic servo- 
mechanism mounted 3.05 m above the wavemaker hinge. The forced motion 
of the wavemaker may be either random or periodic and is activated by either 
an electronic function generator or a digital time sequence synthesized on a 
PDP 11 E10 digital minicomputer through digital-to-analog converters (DAC). 
A description of a unique inverse finite Fourier transform algorithm (FFT) 
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FIG. 1.—Wave Flume Geometry of OSU-WRF 


developed at the OSU-WRF to generate both periodic and random waves of 
the periodic-random type through a minicomputer is given by the first writer 
and Borgman (5). The periodic motions of the hinged wavemaker which were 
used to verify the linear theory wavemaker design curves for periodic waves 
were also synthesized from a stacked FFT digital computer algorithm using 
the theoretical wavemaker gain function, S/H. A more detailed description 
of the wave generating characteristics of the OSU-WRF is given by Sollitt 
and Huber (9). 

The wave heights and periods selected for experimental verification were 
determined from the dimensionless relative wave height-relative water depth 
dissection plane given by Dean (Ref. 3, Fig. 23, p. 39). Wave heights which 
were 25% of the theoretical breaking wave height were selected for verifying 
the linear wave theory design curves. The desired wave height, H, wave period, 
T, and wave flume water depth in the test section, h-A, were input to the 
minicomputer stacked FFT algorithm; and a discrete sinusoidal time sequence 


4 | 
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for the displacement of the wavemaker was synthesized and output to the hydraulic 
servomechanism which activates the hinged wavemaker through the minicomputer 
DAC. This digital simulation procedure also permits the computation of an 
additional complex-valued wavemaker gain function which is the result of the 
servo-mechanism in the hydraulic system by correlating the measured wavemaker 
motion, S,,, with the minicomputer DAC signal for the theoretical wavemaker 
stroke, S ,. The design wave parameters specified for the experimental verification 
are identified in Fig. 2 on the dimensionless relative wave height-relative water 
depth dissection given by Dean (3). 


Relative Depth 


FIG. 2.—Relative Wave Height-Relative Water Depth Dissection Plane Used to Deter- 
mine Theoretical 25% Breaking Wave Heights for Experimental Verification in 
OSU-WRF [from Dean (3)] 


Experimental verification of the linear theory wavemaker design curves in 
the OSU-WRF consisted of verifying the following wavemaker variables: (1) 
The dimensionless gain function, S / H ; (2) the magnitude of the total dimensionless 
hydrodynamic pressure moment on the wavemaker, M’, and the relative phase 
angle between the hydrodynamic pressure moment and the displacement of 
the wavemaker, B,,; or equivalently (3) the propagating, M), and evanescent, 
M!, components of the total dimensionless hydrodynamic pressure moment. 
The dimensionless hydrodynamic pressure moment was selected for verification 


10 
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instead of the dimensionless hydrodynamic pressure force because verification 
of the hydrodynamic pressure moment avoids having to estimate the moment 
arm of the total hydrodynamic pressure force theoretically [compare Gilbert, 
et. al (4)]. 


Experimentat Data Recorpinc 


An example of a visual analog record of the measured data recorded in the 
OSU-WRF water depth at the wavemaker of h = 3.96 m is reproduced in 


! 
| 


FIG. 3.—Sample Visual Analog Data Record (H = 0.55 m, T = 3.45 sec, and h = 
4.42 m) 


Fig. 3. The dimensionless wavemaker gain function was computed from a ratio 
of the measured wavemaker stroke, S,,, to the measured wave height, H. The 
measured wavemaker stroke was recorded through an LVDT mounted just above 
the wavemaker hinge in the dry well behind the wavemaker (vide, Fig. 1). 
This LVDT signal was calibrated by rotating the wavemaker through a discrete 
set of arcs measured at the height of the wavemaker piston and correlating 
the LVDT signal with those measured rotations. The wave heights were computed 
from a measurement of the instantaneous water surface elevation in the test 
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section of the wave flume by a Sonic System Model 86 sonic wave profiler. 
The sonic wave profiler was calibrated by displacing the sonic transducer through 
a discrete sequence of vertical distances measured above the still water level 
in the absence of waves and correlating the sonic transducer signal with these 
measured displacements. 

All of the experimental amplitudes and relative phase angles were computed 
from FFT coefficients which were computed from the digitized analog signals 
using the integer FFT algorithm contained in the SPARTA subroutine of the 
Lab Peripheral System written by DEC for the PDP 11 E10. Exactly four periods 
of each measured signal were analyzed by the FFT algorithm at a sampling 
interval determined from the specified design wave period according to 


FIG. 4.—Schematic Representation of Hinged Wavemaker Piston and Flap System 
in OSU-WRF 


were recorded just after steady-state conditions were observed to be established 
on the visual analog records but before any major reflections could propagate 
back from the terminating beach. The measured data were, therefore, not 
corrected for reflections from the terminating beach. 

The measured analog signals shown in Fig. 3 were Fourier analyzed by FFT 
algorithm in order to more accurately estimate their magnitude and phase [compare 
to Clough and Penzien (2)]. For a complex-valued discrete time sequence {U,}, 
the discrete FFT pair is given by 


(1) 
1,024 
© in which T = the specified design wave period. The data used in the analyses 
eit) 
Li 
Rit) 
wy 
Ped 
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m=0 
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a= n€Xp—i > MHYV, 1, 2, 
Ww 


n=0 


in which the real and imaginary components of the complex-valued FFT 
coefficients are identified by 


X,, = Re{X,,} + ilm{X,,} 


Only the real part of the discrete complex-valued time sequence {U,,} represents 
the physical process synthesized from Eq. 2a for an unstacked FFT algorithm. 
In order for the inverse transform to be exact, we must also require that 


Since exactly four periods of measured data were digitized, the amplitude 
and phase for the fundamental of the periodic signals will be displayed by 
the fourth harmonic of the FFT coefficients. To illustrate, substitute Eq. 1 
into Eq. 4 and obtain for W = 4,096 


Fourier analyzing exactly four periods of the measured data served to average 
the amplitudes. The accuracy of the integer FFT algorithm which partitions 
each analog signal into unsigned integer values between 0-4,095 was verified 
by manual computations using values computed from the calibrated visual analog 
records (e.g., Fig. 3). 

A check of the periodicity in each analog signal was made by comparing 
the magnitude of twice the square modulus of the FFT coefficient for the fourth 
harmonic to the total energy content of the amplitude spectrum for the discrete 
measured time sequence; i.e.: 


in which m,= > |X,,|” 


m=0 


and j = subscript symbol which identifies a particular measured data signal. 
The closer the value of the ratio computed from Eq. 6 is to 100%, the more 
linear or sinusoidal the process. The lower the value of ,, the less the variance 
of the Fourier analyzed signal may be explained by the fourth harmonic. 

One way in which the theoretical total pressure force on a wavemaker could 
be verified experimentally would be by measuring with pressure transducers 
the hydrodynamic pressure at discrete vertical elevations along the hinged 
wavemaker face and then integrating these pressures numerically. On the other 
hand, it is the total force required to drive the wavemaker that is of interest 
in the design of wavemakers; and it was the moment induced by this force 
which was measured to verify the linear wavemaker theory in the OSU-WRF. 

A differential pressure transducer was installed in the MTS 406 Controller 
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System which activates the hinged wavemaker in order to measure the time 
varying pressure signal of the hydraulic oil piston pressure which oscillates 
the wavemaker piston and which is used in the feedback control system (vide, 
Fig. 4). The voltage output from this differential pressure transducer was calibrated 
to provide an analog signal of the hydraulic oil piston pressure force which 
drives the wavemaker piston and oscillates the hinged wavemaker. This measured 
hydraulic oil piston pressure force could then be used with Newton’s second 
law of motion to measure the hydrodynamic wave pressure moment on the 
wavemaker. Again, the hydrodynamic pressure moment was selected for verifica- 
tion of the design curves presented in the companion paper since the moment 
arm to the piston pressure force, H, could be accurately measured and we 
avoid having to estimate theoretically the hydrodynamic moment arm to the 
hydrodynamic pressure force, /. 

The voltage signal from the hydraulic oil differential pressure transducer was 
calibrated by correlating the voltage output from the hydraulic oil pressure 
transducer required to maintain the wavemaker in a neutral vertical position 
against a static column of water of known height on one side only of the 
wavemaker. A calibration curve for the hydraulic oil piston pressure force was 
obtained by increasing the static water depths on one side of the wavemaker 
by flooding the wave flume at discrete increments of water depth. No hydraulic 
oil piston pressure force was required to maintain the wavemaker in a neutral 
vertical position until approx | m of static water pressure moment was attained 
in the wave flume. This provided an estimate of the static Coulomb friction 
force provided by the pressurized plastic seals located on the vertical sides 
of the wavemaker. 

The measured hydrodynamic moment on the wavemaker was calculated from 
Newton’s second law of motion using the calibrated differential hydraulic oil 
piston pressure transducer signal from the hydraulic servo-mechanism recorded 
both with water in the wave flume and with the wave flume drained. The dynamic 
equations of motion required to measure the hydrodynamic pressure moment 
on the OSU-WRF hinged wavemaker are derived in the following section. 


Dynamic Equations 


The dynamic equation for the conservation of horizontal momentum for small 
amplitude oscillations of the wavemaker piston section with water in the wave 
flume shown in Fig. 4 may be approximated by (compare to Ref. 7): 


mS(t) + pS(t) = P(t) + N(t) — R(t) 


in which the dot denotes temporal differentiation of the wavemaker piston 
displacement, S(t); m = mass of the piston; » = viscous resistance of the 
piston [Timoshenko and Young (10)]; P,,(¢) = hydraulic oil pressure force of 
the piston with water in the wave flume; and N(t) and R(t) = the dynamic 
nitrogen and connecting rod reaction forces, respectively. The dynamic equation 
for the conservation of angular momentum for small amplitude motions of the 
hinged wavemaker section with water in the wave flume shown in Fig. 4 may 
be approximated by 


16(t) + v = M(t) — HR(t) 
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in which J = moment of inertia of the hinged wavemaker about the hinge; 

viscous resistance of wavemaker flap; M(t) = the hydrodynamic wave 
pressure moment on the hinged wavemaker; and H = moment arm of the 
connecting rod reaction force. For small amplitude motions of the piston, the 
angular displacement of the wavemaker may be approximated by 


tan 6(¢) = 


in which a positive horizontal displacement of the piston induces a negative 
angular rotation of the hinged wavemaker in the right-handed Cartesian coordinate 


TABLE 1.—Summary of FFT Coefficient Ratios, y,, and Measured Phase Angles for 
Depth at Wavemaker, h = 3.96 m 


Ys Yr, Ys YP, 
stroke pressure stroke pressure B,, in By, in S,,/Sr Bay, in 
(water) (water) (drained) | (drained) | degrees| degrees | (water) | degrees 
(2) (3) (4) (5) (6) (7) (8) (9) 


system employed in the companion paper. For the experimental verification 
in the OSU-WRF, S/2H < 0.16 and the approximation given by Eq. 10 will 
have a maximum error = 0.8%. Substituting Eq. 10 into Eq. 9 and equating 
connecting rod reaction forces, R(t), the coupled equations of motion with 
water in the wave flume reduce to 


M 
AS + BS=P,+N-— 


I 
in which A=m+—> 


—S(t) 
| 
(1) 
0.075 100.0 98.0 100.0 89.3 66.5 138.2 1.01 -12.4 
:; 0.094 100.0 97.9 100.0 77.2 71.9 141.0 0.99 -15.5 
0.107 100.0 97.8 100.0 85.9 67.4 139.9 0.97 -15.7 
0.126 99.9 97.3 100.0 84.2 69.9 137.8 0.97 —14.6 
0.138 100.0 97.6 100.0 82.3 69.5 137.7 1.00 —-14.8 
; 0.157 99.9 97.6 100.0 80.3 74.7 136.0 1.02 —18.2 
0.170 99.9 97.2 100.0 79.5 77.0 134.3 0.99 —20.3 
0.189 99.9 97.0 100.0 78.6 77.6 133.3 0.94 —20.0 
0.220 99.9 95.9 100.0 77.4 76.3 133.4 1.00 —17.8 
0.251 99.9 96.4 100.0 75.1 82.9 132.7 0.98 —24.6 
0.283 99.8 95.6 100.0 74.8 78.9 132.1 0.94 —19.7 
7 0.314 99.8 95.2 100.0 74.1 86.2 130.2 1.02 —25.5 F 
| 0.346 99.7 94.4 100.0 74.2 82.9 129.8 0.90 —24.3 
H 0.377 99.7 93.8 100.0 52.7 83.5 127.0 1.03 —23.8 
f 0.408 99.6 88.3 100.0 75:2 90.5 128.2 0.89 —29.9 
: 0.440 99.6 88.3 99.9 73.8 86.2 127.9 0.98 —24.0 
, 0.471 99.5 86.1 100.0 74.9 96.6 127.3 0.91 —33.2 
0.503 99.4 84.8 99.9 75.3 89.3 124.6 0.96 —25.1 
0.534 99.3 85.2 99.9 58.5 101.4 123.3 0.93 —36.2 
0.566 99.1 87.7 99.9 61.3 92.5 124.2 0.88 —25.6 
0.628 94.5 86.1 99.4 58.4 95.0 119.3 0.87 —26.5 
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Similarly, the coupled translational and rotational dynamic equations of motion 
with the wave flume drained may be expressed by 


AS + BS = P,{t) 
in which P, = hydraulic oil pressure of the piston with the wave flume drained. 


TABLE 2.—Summary of FFT Coefficient Ratios, y,, and Measured Phase Angles for 
Depth at Wavemaker, h = 4.42 m 


Yr, Ys 
pressure stroke pressure B,, in By, in S,,/Sr Bay. in 
(water) (drained) | (drained) | degrees| degrees | (water) | degrees 


(3) (4) (5) (6) (7) (8) (9) 


Solving Eqs. 11 and 13 for the hydraulic oil pressure moment induced by 
the piston with water and with the wave flume drained yields 


M(t) 
= N(t) + AS(t) + BS(t) 


P(t) = AS(t) + BS(t) 


For simple harmonic motions of the wavemaker piston, we may define the 
following sinusoids: 


S(t) = ry exp i (wt + a) 
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v 
2 (water) 
(1) (2) 
0.075 99.9 99.2 100.0 90.3 64.1 142.5 0.99 —10.9 
0.094 100.0 98.6 100.0 88.3 70.8 140.5 1.02 —13.4 
0.107 100.0 97.9 100.0 86.4 70.6 138.9 1.00 —15.6 
0.126 99.9 98.9 100.0 84.2 74.5 137.4 0.97 -17.3 
0.138 99.9 98.5 100.0 81.1 72.9 137.3 0.95 —16.3 
0.157 99.9 98.8 100.0 80.2 72.3 135.6 0.96 —14.7 
0.170 99.9 98.7 100.0 66.0 72.2 135.5 1.00 —15.2 
0.189 99.9 98.6 100.0 77.5 77.4 134.2 1.02 —18.5 
0.220 99.8 97.2 100.0 74.3 81.7 131.3 0.95 —21.8 
0.251 99.7 98.4 100.0 72.3 75.5 131.1 0.96 —17.3 
0.283 99.8 98.3 100.0 71.5 80.8 129.8 1.03 —22.5 
0.314 99.3 97.8 100.0 72.1 81.2 129.1 0.90 —23.4 
0.346 99.7 97.4 100.0 71.4 78.9 128.5 1.00 —20.2 
: 0.377 99.4 89.2 100.0 69.7 86.6 127.3 0.97 —27.6 
0.408 99.3 88.9 100.0 71.5 83.6 125.3 0.89 —23.3 
0.440 99.4 96.4 99.9 48.2 86.5 127.1 1.02 —25.2 
0.471 98.4 96.4 99.9 50.9 93.1 125.4 0.91 —29.0 
0.503 99.2 95.2 100.0 72.7 87.2 123.5 0.97 —23.5 
0.534 96.8 91.8 99.9 72.6 97.4 123.8 0.97 —30.8 
‘ 0.566 97.8 86.1 99.8 70.0 92.8 121.2 0.88 —25.6 
0.628 98.4 93.6 100.0 74.2 104.0 119.9 0.91 —33.2 
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TABLE 3.—Relative Errors for Wavemaker Gain Function, «,, for Dimensionless 
h=3.96 m 


asa 
Measured | Theoretical percentage 


(5) (6) (7) 


= P, expi(wt + a,) 
Pt) = P, exp i (wt + a,) 
M(t) = Mexpi (wt + ay) 


N(t) = N expi (wt + a,) 


H 
H(t) = + a,) 


in which only the real part of the complex phase vector is to be taken. Substituting 
Eq. 15c¢ and temporal derivatives of Eq. 15a into Eq. 146 and equating real 
and imaginary parts, the magnitude of the coupled inertia and viscous resistance 
of the wavemaker system are given by the following: 


S/H 
a H/L | H/H, ee 
(1) (2) (3) 

0.075 0.021 0.236 6.25 1.75 1.70 —2.9 
0.094 0.025 0.249 4.95 1.42 1.47 +3.4 
0.107 0.027 0.253 4.24 1.25 1.35 +7.4 
0.126 0.029 0.249 3.30 1.14 1.21 +5.0 
0.138 0.027 0.228 2.63 1.20 1.14 -5.3 
0.157 0.030 0.238 2.22 1.07 1.04 —2.9 
0.170 0.031 0.238 1.94 0.99 0.98 -1.0 
0.189 0.032 0.233 1.58 0.89 0.92 +3.3 
0.220 0.034 0.238 1.20 0.83 0.82 -1.2 
0.251 0.036 0.240 0.93 0.73 0.75 +2.7 
0.283 0.036 0.237 0.70 0.66 0.69 +4.3 
0.314 0.039 0.220 0.57 0.73 0.65 -12.3 
0.346 0.040 0.241 0.45 0.55 0.62 +11.3 
: 0.377 0.039 0.237 0.34 0.62 0.59 -5.1 
0.408 0.038 0.231 0.27 0.52 0.57 +8.8 
0.440 0.038 0.228 0.22 0.57 0.55 —3.6 
: 0.471 0.038 0.230 0.18 0.51 0.53 +3.8 
0.503 0.037 0.218 0.14 0.55 0.52 —5.8 
0.534 0.040 0.234 0.13 0.49 0.51 +3.9 
0.566 0.040 0.232 0.11 0.46 0.50 +8.0 
0.628 0.038 0.227 0.08 0.45 0.49 +8.2 

; Note: € = (theoretical — measured) /theoretical x 100%; and J = number of evanescent 

2P, 
Sw 
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Hydrodynamic Moment, < ,,, and for Relative Phase Angle, « p- for Depth at Wavemaker, 


B,,. in degrees 


€,,asa 
Theoretical | percentage percentage 


(9) (10) (13) 


4 
4 
4 
4 
4 
5 
5 
5 
3 
3 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 


in which the drained phase angle B, = a, — a must be a second quadrant 
phase angle for positive semidefinite values of the system inertia and viscous 
resistance [cf. Webster (12) for a review of the phase lag for a damped harmonic 
oscillator with negligible stiffness compared to inertia]. 

Substituting Eqs. 146, 15b, 15d, and 1Se into Eq. 14a and equating real and 
imaginary parts yields the following expressions for the magnitude of the 
hydrodynamic pressure moment on the hinged wavemaker: 


M =H + + N’-2P,P, cos (6, —B,) + 2N[P,, cos (B,, — By) 
— P, cos — By)]}'” 
and for the relative phase angle for the hydrodynamic pressure moment: 
P, sinB, — P, sinB, + N sinBy 
B, = arctan 
P,, cos B,, — P, cosB, + N cos By 


in which B,, = a, — a; and By, =a, —a. 


Mea- 
sured | 
0.4363 0.4077 —7.0 85.8 86.3 +0.6 
0.4036 0.3949 —2.2 84.4 86.3 +2.2 
0.3618 0.3864 +6.4 86.8 86.4 —0.5 
0.3769 0.3737 87.4 86.6 -1.0 
0.3563 0.3653 +2.5 88.2 86.8 -1.7 
0.3513 0.3527 +0.4 88.9 87.1 —2.1 
0.3251 0.3445 +5.6 86.0 87.4 +1.6 
0.3315 0.3322 +0.2 87.4 87.8 +0.5 
0.3051 0.3123 +2.3 87.6 88.5 +1.0 
0.2949 0.2932 —0.6 85.3 89.0 +4.1 
0.2633 0.2753 +4.4 87.6 89.2 +1.9 
0.2469 0.2586 +4.6 84.3 89.1 +5.4 
0.2414 0.2434 +0.8 89.7 88.6 —1.2 eS 
0.2327 0.2296 -1.3 89.7 87.5 —2.5 
0.2358 0.2180 -8.1 82.2 86.0 +4.4 
0.2238 0.2071 -8.1 89.1 84.0 -6.1 
0.2269 0.1978 —14.7 75.3 81.5 +7.7 
0.2193 0.1905 -15.1 88.0 78.6 -11.9 
0.2092 0.1842 —13.6 67.6 75.4 +10.4 
0.2054 0.1795 —14.5 85.7 71.8 —19.3 
0.2087 0.1746 —19.5 81.9 64.0 —27.8 
required for 0.1% accuracy. : 
Sw 
(17a) 
(17b) 
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TABLE 4.—Relative Errors for Wavemaker Gain Function, «., for Dimensionless 
h=442m 


asa 
Measured | Theoretical percentage 


(5) 


Note: € = (theoretical — measured) /theoretical x 100%; and J = number of evanescent 


In the OSU-WRF, the time varying signals from the hydraulic oil piston 
pressures, P(t) and P,(t), along with the wavemaker piston displacement, S(t), 
were recorded and were Fourier analyzed by an FFT algorithm to compute 
the magnitudes and phases for these two signals which are required in Eqs. 
17(a) and 17(b). The time dependent nitrogen gas pressure term, N(t), was 
not recorded directly but had to be estimated from the measured piston 
displacement as described in the following. 

The nitrogen gas system in the MTS designed OSU-WRF is used to backpressure 
the wavemaker piston in order to compensate for the static column of fluid 
on the wave flume side of the hinged wavemaker (vide Fig. 4). The time varying 
nitrogen gas pressure is a result of the small volume changes induced by the 
wavemaker piston stroke. The nitrogen gas pressure on the wavemaker piston 
may, therefore, be approximated by a static plus a time varying component 


according to the real gas equation: 
nZRTA y 

N=N,+N(t)=——_ 
Vy +AyS(t) 


in which n = mass of the nitrogen gas; Z = the compressibility coefficient 
for a real gas; R = gas constant for nitrogen; T = temperature; V, = static 


S/H 
see 
2 H/L | H/H, 

(1) (2) (3) a ia (7) 

0.075 0.023 0.257 6.76 1.44 1.50 4.0 

0.094 0.024 0.244 4.85 1.35 1.35 0.0 

0.107 0.025 0.233 3.91 1.29 1.24 —4.0 

0.126 0.028 0.246 3.26 1.05 1.11 5.4 

0.138 0.030 0.256 2.92 0.92 1.04 1.5 

0.157 0.031 0.243 2.27 0.91 0.95 4.2 

0.170 0.031 0.240 1.97 0.90 0.90 0.0 

0.188 0.031 0.228 1.55 0.90 0.84 —7.1 

0.220 0.035 0.242 1.22 0.70 0.75 6.7 

0.251 0.035 0.234 0.90 0.67 0.68 1.5 

0.283 0.037 0.242 0.71 0.64 0.63 —1.6 

f 0.314 0.038 0.242 0.56 0.53 0.59 10.2 

0.346 0.037 0.231 0.42 0.58 0.56 —3.6 

0.377 0.038 0.234 0.34 0.53 0.53 0.0 

0.408 0.040 0.242 0.28 0.45 0.51 11.8 

0.440 0.039 0.231 0.22 0.52 0.50 —4.0 

; 0.471 0.036 0.218 0.17 0.48 0.48 0.0 

0.503 0.038 0.230 0.15 0.47 0.47 0.0 

0.534 0.039 0.233 0.13 0.46 0.46 0.0 

0.565 0.037 0.220 0.10 0.43 0.45 4.4 

0.628 0.037 0.219 0.07 0.43 0.44 2.3 
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Hydrodynamic Moment, < ,,, and for Relative Phase Angle, <,, for Depth at Wavemaker, 


B,,. in degrees 


Ey, asa €,, asa 
Theoretical | percentage Theoretical | percentage 
(9) (10) (12) (13) 


4 
4 
4 
4 
4 
4 
5 
3 
3 
3 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 


terms required for 0.1% accuracy. 


volume of nitrogen gas; and A, = surface area of the terminal end of the 
wavemaker ram normal to the wavemaker stroke (which is not equal to the 
hydraulic oil piston surface area). 

For the experimental verifications in the OSU-WRF, the ratio of A , S(t)/ Vy 
= 0.16; and we may expand Eq. 18 by the binomial theorem to obtain, 
approximately 


nZRTA y Ay S(t) (A050 ) 


N N N 


in which the static component exactly balances the hydrostatic moment on the 
wavemaker and the dynamic component acts like a stiffness element in opposing 
positive displacements of the piston. The difficulty in recording the dynamic 
nitrogen pressure force directly may be seen from the ratio of the dynamic 
to static component; i.e.: 

N, Vy 


Pressure instrumentation with adequate resolution for dynamic pressures which 


M’ 
Mea- 
(8) i 
0.4191 0.4116 —1.8 84.0 86.3 +2.7 
0.4067 0.3987 —2.0 76.6 86.3 +11.3 
0.4212 0.3902 -7.9 80.8 86.4 +6.4 
0.3976 0.3774 —5.4 81.5 86.6 +5.9 
0.3364 0.3690 +8.8 83.3 86.7 +4.0 
0.3576 0.3565 —0.3 87.0 87.1 +0.1 
0.3356 0.3482 +3.6 85.6 87.3 +2.0 
0.3682 0.3360 —9.6 81.4 87.7 +7.2 
0.3138 0.3162 +0.8 79.4 88.4 +10.2 
0.2967 0.2973 +0.2 88.0 88.9 +1.0 
0.2804 0.2795 —0.3 81.0 89.2 +9.2 
0.2596 0.2629 +1.3 85.7 89.2 +4.0 
0.2522 0.2477 —1.8 86.1 88.8 +3.0 
0.2458 0.2340 —5.0 78.3 87.9 +11.0 . 
. 0.2138 0.2224 +3.9 84.6 86.6 +2.4 
0.2183 0.2112 —3.4 78.1 84.9 +8.0 
0.2303 0.2016 —14.3 73.5 82.6 +111 
0.1965 0.1939 -1.3 78.9 80.0 +1.4 
0.2005 0.1871 —7.2 65.5 77.0 +14.9 
0.1947 0.1818 —7.1 74.1 73.6 —0.5 
0.1982 0.1756 —12.9 56.8 66.3 +14.4 
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are very much less than 16% of the static component are difficult to obtain 
economically. Consequently, the experimental verifications were analyzed by 
replacing the exact dynamic nitrogen gas terms in Eqs. 17(a) and 17(b) with 
an equivalent linear spring element approximated by 


in which the equivalent spring constant, K, is determined by a static moment 
balance between the hydrostatic pressure moment on the wavemaker and the 
static nitrogen gas pressure according to 


N 
BA ,(h —8)° 
N 


in which B = width of the hinged wavemaker. Substituting Eq. 21 into Eqs. 


FIG. 5.—Booker (1) Phase Vector Diagram for Dynamic Oscillations of OSU-WRF 
Hinged Wavemaker System 


17(a) and 17(b) and evaluating the equivalent nitrogen gas spring constant, K, 
from the static water depth in the wave flume according to Eq. 22, the following 
relationships for the magnitude and phase angle, respectively, for the wavemaker 
hydrodynamic pressure moment are obtained: 


2 
M=H + P2+ («=) — 2P,,P, cos ®, —B,,) — KS[P, cosB,, 
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P,, sin B,, — P,sinB , 
6 = arctan 


P,, cos B,, — P, cos B, —- K — 


since =ay-a=a—a=0. 
Experimenta Resutts 


Experimental measures for the wavemaker gain function, S/H; magnitude 
of the dimensionless hydrodynamic pressure moment, M’; and relative phase 
angle, B,,, were obtained for 21 values of relative wave frequency for each 


hyn) sym; 
| 


« | 


FIG. 6—Comparison of Measured Wavemaker Gain Function, S_, / H, with Theoretical 
Wavemaker Gain Function, S,/H, for Depths at Wavemaker, h = 3.96 m and = 
442m 


of two different water depths which cover one decade of relative water depth, 
(h — A)/L,. Tables 1 and 2 summarize the dimensionless FFT coefficient ratios, 
,, computed from Eq. 6 for the experimental data recorded in the two water 
depths. These dimensionless FFT coefficient ratios indicate that the data measured 
for the wavemaker stroke with water (Tables 1 and 2, Cols. 2) and with the 
wave flume drained (Tables 1 and 2, Cols. 4) are very well represented by 
a single linear sinusoid with an amplitude given by the magnitude of the FFT 
coefficient for the fourth harmonic. The hydraulic oil piston pressure force 
ratio with water in the wave flume which are summarized in Tables 1 and 
2, Cols. 3, are also well explained by the magnitude of the FFT coefficient 
for the fourth harmonic except for dimensionless relative frequencies, 2, greater 
than 0.314 which are deep-water wave conditions. On the other hand, the 
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magnitudes for the hydraulic oil piston pressure force with the wave flume 
drained which are summarized in Tables | and 2, Cols. 5, are not as well 
correlated with the magnitudes of the FFT coefficients for the fourth harmonic 
and this lower correlation continues to decrease with increasing relative wave 
frequency, . This explains why the relative phase angles with the wave flume 
drained, B,, tabulated in Tables 1 and 2, Cols. 7, rotate toward 90° instead 
of 180° with increasing values of dimensionless wave frequency, 2, as would 
be expected from an inertially dominated dynamic system with linear viscous 
resistance [compare to Timoshenko and Young (10) or Webster (12)] . The reasons 
for this counter rotation may be observed from the hydraulic oil piston pressure 
signal, P,, shown on the sample analog trace in Fig. 3. The high frequency 
oscillation which may be seen superimposed on the fundamental hydraulic oil 
pressure signal is a result of the feedback control system and its amplitude 
is a greater percentage of the fundamental amplitude for, P,, than it is for 


a. 
r= 
=| 
@ 
o 
Q 


= 0.7692 
8/h= 0.0192 
A/sh= 0.2308 


FIG. 7.—Comparison of Measured and Theoretical Dimensionless Hydrodynamic 
Wavemaker Amplitudes, M’ and M/, for Depth at Wavemaker, h = 3.96 m 


P,,. This high frequency error in the drained piston oil pressure force precludes 
an accurate estimate of the system inertia and viscous resistance via Eqs. 16(a) 
and 16(5); however, it does permit an approximate evaluation of the hydrodynamic 
pressure moment on the wavemaker. 

Cols. 8 and 9 of Tables 1 and 2 summarize the magnitudes of the wavemaker 
relative gain function, S,,/S ,, and relative phase angles, B,,, which were computed 
by correlating the measured wavemaker stroke from the LVDT signal, S,,, 
with the theoretical wavemaker stroke synthesized from the minicomputer FFT 
algorithm, S,. The increase in the magnitude of the phase lag of the measured 
wavemaker stroke behind the theoretical wavemaker stroke and the decrease 
in the magnitude ratio may be observed from Tables | and 2, Cols. 9 and 
8 to depend upon the dimensionless relative wave frequency, 2. 

Tables 3 and 4 summarize the dimensionless magnitudes and relative phase 
angles for the hydrodynamic pressure moments measured on the OSU-WRF 
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hinged wavemaker which were computed from the measured FFT coefficients 
by Eqs. 23(a) and (b). Relative measures of the degree of nonlinearities in 
the measured wave data are provided by: (1) The relative wave steepness values, 
H/L, tabulated in Col. 2 using the measured wave height and linear wave 
theory estimates for the wave length; (2) the ratio of measured wave height, 
H, to the theoretical breaking wave height, H,, tabulated in Col. 3; and (3) 


the Ursell parameter, U [Ursell (11)] , tabulated in Col. 4, which may be computed 
from the following ratios: 


Dimensionless relative errors are tabulated for each of the three dimensionless 
wavemaker variables measured from 


( theoretical — measured ) 
theoretical 


Negative values for these relative errors indicate an underprediction by linear 
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FIG. 8.—Comparison of Measured and Theoretical Dimensionless Hydrodynamic 
Wavemaker Amplitudes, M’ and M/, for Depth at Wavemaker, h = 4.42 m 


wave theory derived in the companion paper while positive values of relative 
errors indicate overprediction by linear wave theory. Col. 14 lists the number 
of evanescent eigenmode terms required in the infinite summation in order to 
obtain an accuracy of 0.1%. 

Fig. 5 is a schematic Booker phase vector diagram (1) for a typical set of 
dynamic oscillation signals recorded for the OSU-WRF hinged wavemaker system. 
This type of diagram is quite powerful for quickly checking the relative accuracy 
of approximate analytical dynamic structural models. An a priori analyses of 
the measured data in the OSU-WRF assumed that the ratio of the magnitude 
of the dynamic nitrogen gas pressure force to the magnitude of the drained 
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piston oil pressure force, N/P,, would be negligible and was omitted from 
initial computations. A quick series of Booker phase vector diagrams using 
the FFT magnitude and phases for P,, and P, quickly demonstrated that the 
magnitudes of the resultant hydrodynamic moment vector were consistently 
much larger than values estimated by linear wave theory and that the relative 
phase angles were consistently much smaller than values predicted by linear 
wave theory (vide Fig. 5 and omit the phase vector—KS). It was readily visualized 
as a result of the Booker phase vector diagrams that the vector addition of 
a dynamic phase vector which was in-phase with the inertia of the wavemaker 
system (i.e., in the negative reference direction in Fig. 5) would yield a resultant 
phase vector for the hydrodynamic pressure moment with approximate magnitude 
and phase values closer to those predicted by linear wave theory. The equivalent 
linear spring model for the dynamic nitrogen gas pressure force was then 
reintroduced into the analyses and a posterior check of the ratio KS /P, verified 


120 


Hin 
0.7692 0.0192 0.2308 
© 0.6897 0.0172 0.2069 


FIG. 9.—Comparison of Measured and Theoretical Relative Phase Angles for Wave- 
maker Hydrodynamic Pressure Moment for Depths at Wavemaker, A = 3.96 m (—— 
Solid) and 4.42 m (---- Dashed) 


that the ratio could not be neglected in the analyses. The Booker phase vector 
diagram was quite powerful for rapidly analyzing multiple oscillatory signals 
and for evaluating approximate analytical dynamic structural models. 

Fig. 6 compares the measured wavemaker gain function values, S,,/H, with 
the theoretical values computed in the companion paper. The absolute values 
of the maximum relative errors are less than 12.5% for the 42 data values 
measured. These values for the relative errors for the dimensionless wavemaker 
gain function are consistent with values previously reported by others for hinged 
wavemakers. 

Figs. 7 and 8 compare measured and theoretical values for the magnitudes 


of the dimensionless wavemaker pressure moment for the propagating mode 
which were computed from 


= M’ sinB, 
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and for the evanescent mode which were computed from 
M! = M’ cosBy 


The absolute value for the maximum relative error, €,,, is approx 20% (Table 
3, Col. 10). The magnitudes of the relative errors for the dimensionless hydro- 
dynamic pressure moment on the hinged wavemaker increase with increasing 
relative wave frequency, 2. These measured magnitudes of wave pressure moment 
on a hinged wavemaker appear to be unique. 

Fig. 9 graphically compares the measured and theoretical relative phase angles, 
B,,, for the hydrodynamic pressure moment on the hinged wavemaker. The 
magnitudes of the relative errors for the hydrodynamic moment phase angle, 
B, are less than 28% (Table 3, Col. 13) and, again, increase in magnitude 
with increasing relative wave frequency, 2. 


Conctusions 


The dimensionless design curves for hinged wavemakers of variable draft 
for wave flumes which consist of two constant depth domains separated by 
a gradually sloping transition region are compared with unique experimental 
data recorded in the OSU-WRF in each of two water depths for 21 values 
of dimensionless wave frequency and relative wave heights which are approx 
25% of the theoretical breaking limit. The magnitude and phase angle of the 
hydrodynamic pressure moment on the wavemaker were estimated from the 
coupled translational and rotational dynamic equations of motion by FFT 
coefficients which were computed from the wavemaker stroke and from the 
measured hydraulic oil piston pressure force with water and with the wave 
flume drained. The FFT coefficients for the fourth harmonic demonstrated good 
correlation for a linear process except for the drained wave flume condition 
in which a high frequency component in the piston pressure force which was 
induced by the feedback control system diminished the percent contribution 
by the fundamental periodic component to the total energy content of the amplitude 
spectrum. The dynamic component of the nitrogen gas pressure force was 
approximated by an equivalent linear spring in which the linear spring constant 
was evaluated from the real gas equation and a static moment balance between 
the nitrogen gas backpressure and the static column of fluid in the wave flume. 
The Booker phase vector diagram (1) proved to be quite powerful for rapidly 
comparing measured FFT phase vectors of oscillatory signals with approximate 
analytical dynamic structural models in posterior analyses. 

While measured data for the dimensionless wavemaker gain function, S/H, 
have been previously published, the comparisons with the magnitude and phase 
angle of the hydrodynamic pressure moment appear to be unique. 
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Aprennix Il.—Nortation 
The following symbols are used in this paper: 


coupled masses of piston and wavemaker; 


Ay = surface area of nitrogen piston; 
B = total width of wavemaker; 7 
B = coupled mechanical damping of piston and wave- 


maker; 

gravitational constant; 

deterministic wave height; 

still water depth at wavemaker; 

moment of inertia of hinged wavemaker about hinge; 
imaginary unit; 

total number of evanescent eigenmodes required to 
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obtain 1% accuracy in infinite summations; 

wave number for propagating eigenmode in constant 
depth domain, h, adjacent to wavemaker; 

wave number for evanescent eigenmodes (n = 2); 
linear theory wave length; 

deep water wave length; 

moment arm to hydrodynamic pressure force; 

time dependent (amplitude) dimensionless total hydro- 
dynamic moment on the wavemaker; 

mass of hydraulic piston which oscillates hinged 
wavemaker; 

zeroeth spectral moment (= variance of time se- 
quence); 

time dependent (amplitude) nitrogen gas pressure 
force; 

mass of nitrogen gas; 

summation indices; 

time dependent (amplitude) hydraulic oil piston pres- 
sure force with wave flume drained; 

time dependent (amplitude) hydraulic oil piston pres- 
sure force with water in wave flume; 

time dependent connecting rod reaction force; 

gas constant; 

time dependent (double amplitude) displacement of 
wavemaker stroke measured at height of wavemaker 
piston; 

wave period; 

temperature of nitrogen gas; 

time; 

complex-valued discrete time sequence; 

Ursell parameter; 

volume of nitrogen gas required to back-pressure 
piston; 

integer number of values in discrete time sequence; 
complex-valued FFT coefficient; 

horizontal coordinate axis with origin located at undis- 
turbed water level at wavemaker; 

vertical coordinate axis with origin located at still water 
level at wavemaker; 

compressibility coefficient for real gas; 

arbitrary initial phase angle for wavemaker variable, 
J; 

relative phase angle measured between wavemaker 
variable, j, and displacement of wavemaker; 

ratio of twice square modulus of FFT coefficient for 
fourth harmonic for wavemaker variable, j, to total 
energy content of amplitude spectrum; 

vertical rise in depth between two constant depth wave 
flume domains; 
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equal discrete intervals of frequency used in FFT 
digitizing process; 

equal discrete intervals of time used in FFT digitizing 
process; 


height of wavemaker hinge measured above bottom; 


dimensionless relative error; 

height of wavemaker piston measured vertically above 
wavemaker hinge; 

instantaneous water surface elevation measured posi- 
tive upwards from undisturbed still water level; 
angular rotation of wavemaker about hinge; 
equivalent spring stiffness for dynamic nitrogen gas 
force; 

wave number for propagating eigenmode in constant 
depth test domain, h — A; 

viscous resistance of piston; 

viscous resistance of wavemaker; 

fluid density; 

relative water depth; and 
radial frequency (= 27/T). 


dimensionless variables; and 
temporal derivative. 


theoretical breaking wave height; 

with wave flume drained; 

evanescent eigenmode (n = 2); 

wavemaker gain function; 

hydrodynamic pressure moment; 

measured value; 

nitrogen gas pressure force; 

deep water wave conditions; 

propagating eigenmode; 

wavemaker stroke; 

theoretical wavemaker stroke (discrete time sequence 
synthesized by minicomputer algorithm); and 

with water in wave flume. 
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NUMERICAL ANALYSIS OF FLOW 
IN SEDIMENTATION BASINS 


By David R. Schamber,’ A. M. ASCE and Bruce E. Larock,” M. ASCE 


INTRODUCTION 


Primary sedimentation by gravity in large basins and tanks has been an integral 
part of most major water and sewage treatment plants for many decades. Current 
understanding of the process is incomplete, however, and is impeded primarily 
by a lack of detailed knowledge of the velocity field in the basin. In partial 
response to this situation, this article describes a numerical model which predicts 
the velocity field in sedimentation basins. Since the flow is turbulent within 
these tanks, determination of the velocity field is difficult. The point velocities 
within the basin are determined by the finite element solution of five coupled, 
nonlinear partial differential equations. The structure of the turbulence is 
represented by an effective or ‘“‘eddy’’ viscosity model which depends on the 
turbulent kinetic energy and its rate of dissipation. Once the velocity field is 
determined, the distribution of particle concentration throughout the basin could 
then be obtained by solution of a linear convection-diffusion equation. 


BackGrounb 


Since primary sedimentation by gravity often precedes other unit operations, 
the performance of all processes could potentially be influenced by the efficiency 
of the settling basin. The writers view the design of basins for increased particle 
removal efficiency as primarily a fluid mechanics problem since the basin velocity 
field strongly influences the path taken by particles entering the tank. The ability 
to test basin performance for various tank geometries, different inlet and outlet 
structures, and various loading conditions and thermal regimes is essential for 
determining optimal design criteria. Determining flow patterns in basins is the 
key to evaluating this performance. Flow patterns may be generated by experi- 
mental model studies, both direct field measurement and full-scale experimental 
studies, or mathematical computer modeling. Experimental studies, while useful, 
are expensive and relatively inflexible when the goal is to examine a variety 
of tank geometries and various inlet and outlet structures. Direct field measure- 
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ment is useful for calibration and checking of computer models but by its very 
nature is not predictive. Therefore the writers direct their efforts toward numerical 
solutions of the equations governing the flow in settling basins as the most 
productive long-term alternative. 

The large scale of most tanks [typical vertical and horizontal dimensions 
of 10 ft (3 m) and 100 ft (30 m) respectively], together with the small kinematic 
viscosity of water [v ~ .00001 sq ft/sec (.0000009 m7’/s)], render the flow 
turbulent even though the velocities are in the range of 1 ft/min-10 ft/min 
(0.3 m/min-3.0 m/min). With these characteristic values the tank Reynolds 
number (depth x velocity/v) ranges from 17,000-170,000. These values agree 
well with values computed by Larsen (14) from experimental measurements 
in secondary clarifiers. Larsen also found turbulence intensities to range from 
10%-20% (relative to the maximum mean value of velocities in the bottom 
current) throughout the basin, although even higher values were found in the 
inlet zone. 

The general circulation pattern in rectangular and circular primary sedimentation 
basins, including the effects of surface wind, variable influent temperature, 
solar heating, and unsteady loading conditions, is time variant (unsteady) and 
three dimensional. The cost of a transient three-dimensional computation is 
currently prohibitively expensive in terms of computer execution time and storage. 
However, Montens (20) has shown by direct measurement on full-scale circular 
basins that the flow is very nearly radial in the absence of surface wind. In 
rectangular basins the flow remains three-dimensional due to corner and side 
effects; however, the flow in a vertical plane taken on the tank center line 
from inlet to outlet may reasonably be expected to be nearly two dimensional 
in the absence of wind stresses. 

Thermal effects may induce density differences within the tank and cause 
the flow pattern to stratify. Although the inclusion of temperature effects into 
the computational algorithm appears relatively straightforward, it also increases 
by at least three the number of partial differential equations which must be 
solved (24). The proposed model is therefore limited in application to those 
Operating conditions where the influent stream and fluid temperature within 
the tank are nearly equivalent. 

The sediment concentration in the influent is assumed not to affect the fluid 
density because its magnitude is typically .0002 by weight (19). Presumably, 
this assumption would be invalid for certain basin operating conditions. The 
flow pattern established by an influent stream with little fluid inertia (correspond- 
ing to low loading rates) would be influenced by small concentration differences. 
An inlet stream with density greater than that of the fluid within the basin 
would plunge to the basin floor. By contrast, an inlet jet of lighter density 
would move as a layer over the top of the tank. Low inflow velocities would 
also produce weak turbulence and the high turbulence model (applicable for 
high loading rates) adopted for this study must then be modified appropriately 
(11). The proposed model is therefore restricted to those flow conditions which 
allow the hydrodynamic model to be uncoupled from the transport equation 
for particle concentration and the thermal energy equation. Primary sedimentation 
by gravity of a dilute suspension of discrete particles should adequately be 
represented by such a hydrodynamic model. However, the proper analysis of 
secondary clarifier performance, where density differences are significant, would 
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require the mathematical description of stratification caused by species con- 
centration. Hamlin and Wahab (9) consider further the effects of temperature 
and concentration differences. 

The scope of this paper is therefore limited to two-dimensional steady, 
nonstratified flow in the x-z (horizontal-vertical) or r-z (radial-vertical) plane. 


Matuematicat Mover 


The equations describing two-dimensional, steady, turbulent, nonstratified flow 
in a rectangular or circular sedimentation basin are (25): 


—(r"U) +— (r" W) = 0 
or dz 


aU oU om a 1 
‘w)+—(P+ uu) + — (uw) =0 (2) 
or or az Fo ar 


ow W 42 a 
U—— + W— +—(uw) +— (P+ ww) =0 

or 0z or 0z 
Equation 1 expresses conservation of mass, while Eqs. 2 and 3 express 
conservation of momentum in the horizontal/radial and vertical directions, 


Free surface Over flow 


Solid | 
boundaries— 


FIG. 1.—Schematic Tank Configuration 


respectively. For m = 1 these equations describe axisymmetric flow in the 
r-z plane, and for m = 0 they describe planar flow in the x-z Cartesian plane. 
The coordinate r represents the x coordinate when m = 0. Variables in Eqs. 
1-3 have been nondimensionalized on a characteristic length A, and characteristic 
velocity U,, and F? = U2/(gh,) is a representative Froude number with g = 
ratio of weight to mass. The independent and dependent variables are defined 
with reference to Fig. 1 as follows (starred and characteristic variables are 
pommel r* = rh, = horizontal coordinate; z* = zh, = vertical coordinate; 
U* = UU, = time-averaged mean horizontal or radial velocity; w* = WU, 

= time- -averaged mean vertical velocity; u*u* = uu U2, w*w* = ww U3, 
w* = uw U2 are Reynolds stress terms. The overbar denotes time averaging, 
and u and w are instantaneous velocity fluctuations about their respective means. 
The nondimensional pressure P = P*/(pU3) in which p = fluid density, is 
the deviation from the nondimensional hydrostatic pressure, i.e., the actual 
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pressure is given by P + (h — z)/F2 in which h = h*/h, is the free surface 
elevation. The viscous terms which normally appear in Eqs. 2 and 3 are omitted 
here because turbulent transport effectively swamps molecular transport. Inspec- 
tion of Eqs. 1-3 reveals that the number of unknowns exceeds the number 
of equations—this is the familar closure problem associated with turbulent flows. 
Numerous closure models of varying complexity have been proposed in the 
literature (4,16,17). An outline of several models with finite element solutions 
to some one-dimensional flow examples is presented in Ref. 26. 

The kinetic energy-dissipation or k-e turbulence model which determines the 
variation of the turbulent eddy viscosity throughout the basin from the solution 
of two turbulence transport equations will be used to close the present mathemati- 
cal model. This turbulence model has been tested and experimentally verified 
for a number of free-shear flows, boundary layer flows, internal flows, and 
flows with recirculation (15,17). More recently the model has been successfully 
used to model several complex recirculating flows (8,12,21). These computational 
results also compare favorably with available experimental data. It is therefore 
expected that the k-e model will also perform adequately for sedimentation 
basin computations where several of these same flow features coexist. 

By direct analogy to laminar flow, Boussinesq suggested that the turbulent — 
stresses could be represented by expressions involving products of mean velocity 


gradients and the turbulent kinematic viscosity v,. This leads to the isotropic 
closure 


in which the turbulence kinetic energy k = (1/2)(uu + vv + ww). Evidently 
v, must be several orders of magnitude larger than the molecular kinematic 
viscosity v to conform to the observed ability of the turbulent motion to enhance 
momentum transport. Unlike the molecular viscosity, v, is not a property of 
the fluid but varies from point to point within the flow field, reflecting the 
local structure of turbulence. 

Computations of turbulent flow with v, = constant (or known specified variation) 
have been relatively successful (13) when the magnitude of v, is known a priori 
and have the practical advantage that no additional differential equations are 
required for this turbulence closure. Since the k-e model computes the distribution 
of v, from a simulated turbulent flow, one does not require a priori knowledge 
of the distribution of v,. By dimensional considerations v, is proportional to 


the kinetic energy of turbulence and its rate of dissipation ¢€ in the following 
way: 


au 
uu = —k — 2v,— (4) 
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ww =—k — 2v,—— (6) 
| 3 az 
au aw 
0z or 
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Here = v(au,/ax Ou, / ax is the isotropic dissipation rate. The variables 
k and € are computed from transport equations (25) which have the following 
form in steady two-dimensional flow 


ok 
+ 
0z 


The diffusion of k and e, respectively, are given by 


ér \r"o, ar dz \o, 
+- 

ar \r"o, ar az \o, dz 


(7) 


The constants appearing in Eqs. 8, 10, 13, and 14 have the following values 
(17): ¢, = 0.09, c,, = 1.44, c., = 1.92, o, = 1.0, and o, = 1.3. The constants 
Co, ¢,, and o, may be inferred from experimental data (27), while c,, and 
o, are determined by computer optimization, i.e., many calculations were 
performed in which the constant was systematically varied, and the value which 
gave the best overall agreement with experimental results was selected. Variables 
appearing in Eqs. 4-14 are nondimensional such that k* = kU? , «* = (€U2)/(ho), 
and v* = v, U,hy. 

Equations 1-14 represent a closed system of five coupled, nonlinear partial 
differential equations. The five dependent variables are U, W, P, k, and « 
while the independent variables are r and z. Pertinent boundary conditions are 
required to yield a unique solution. These conditions are described later. 


Finite Evement Sovution of Matuematicat Mover 


Elliptic boundary value problems of the type presented herein are generally 
solved by finite difference or finite element techniques. Launder and Spalding 
(17) originally solved the k-e model for a variety of flows by use of the 
Patankar-Spalding (22) finite difference algorithm which is basically a marching 
scheme applied to a parabolic equation set. More recently elliptic problems 
involving recirculation have also been solved with finite differences (12,21). 

In the present work the equations are cast in integral form by the Galerkin 
method of weighted residuals and solved by the finite element method. This 


HY5 579 
k? 
ok 
= or 
€ 
U— + W— —F— 6,2 — es AD 
or az k 
in which the production term F is 


580 MAY 1981 HY5 


approach was chosen because it allows one in principle to use a nonuniform 
computational mesh, to seek solutions in geometrically complex nonpolygonal 
domains, and to avoid difficulties in applying boundary conditions on such domain 
boundaries. An added advantage of this approach, as opposed to upwind or 
penalty function methods, is that the presence of suspicious oscillations in a 
solution appears to be a clear indicator that the computational mesh is locally 
too coarse and should be refined (7). 

Expressions for the equation residuals f, written for an element e having 
domain 2, and a boundary I ,, are (25): 


ow aw 
Sw= | w=) an- | | 
a, or a 


— an 
+(P+ + | [uw 1, + (P+ ww)/,] 
Zz 


ak ak 
A= U— + dn+ 
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ak aN, 
- r 
Oz r, 


de ON, 
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oe 1.) 
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in which the velocity correlations or Reynolds stresses are defined by Eqs. 
4-7. The terms M, and N, are the appropriate basis functions for the dependent 
variables. All second derivative terms in the transport equations have been 
integrated once by parts so that C° elements may be used. The pressure term 
has also been integrated by parts to facilitate the application of normal-stress 
boundary conditions. The line integrals are evaluated only where I’, coincides 
with the global domain boundary I; and /, and /, = direction cosines of the 
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unit outer normal to I’. The term (dh/@r) in Eq. 2 is omitted in Eq. 16 because 
of use of the rigid lid approximation (mentioned with boundary conditions). 

During the 1970s much attention has focused on the correct order of approx- 
imation for the primary variables appearing in the equations. Currently a linear 
variation M, for pressure and a quadratic variation N, for all other variables 
is assumed. While it is now becoming clear why the pressure must be approximated 
one order lower than the velocity components (23), the same cannot yet be 
said with surety for k and e. Eight-node quadrilateral elements are used in 
this work with no known difficulties that can be traced to element choice. 


Numenicat So.ution 


The nonlinear equation set resulting from the finite element discretization 
process is solved by a variation of Newton’s method. Experience with one-dimen- 
sional flow problems (26) indicates that a major difficulty in applying Newton’s 
method is the proper initialization of the solution vector. Another difficulty 
is the sheer size of the equation set when the solution for all five variables 
(U, W, P,k,«€) is sought simultaneously. 

To form the global equation set, one must appropriately sum the elemental 
contributions represented by Eqs. 15-19. The resulting set of nonlinear equations 
may be compactly written as 


f(x) = 0 


in which f = vector representing the system of equations; and x = vector 
of nodal unknowns for the system. 


The basic Newton scheme solves the equation set by writing 
A(x") 5x” = —f(x") 


in which n = iteration counter; A = Jacobian of equation matrix; and 5x” 
= set of incremental corrections to be made to the vector x”. Solution of Equation 
21 gives the following new estimate to the solution 


x 


in which a <= | = a relaxation parameter. 

In this work the mean flow equations (Eqs. 15-17) and the turbulence equations 
(Eqs. 18-19) are solved alternately rather than simultaneously. This method 
is easier to initialize and reduces core storage requirements substantially. 
However, more iterations are required to obtain a converged solution. To start, 
this approach requires a specified initial distribution of k and € (constant values 
or simple linear variations of k and « are generally used) which yield a reasonable 
distribution of v, for the given problem. The mean flow equations are then 
solved by Newton’s method with U, W, and P initialized with constant values. 
For example, P (which is the deviation from the hydrostatic pressure) may 
be set to zero. The U component may be assigned a value equal to the flow 
rate divided by the cross-sectional area, and the W component may be set 
to zero. 

Having generated a solution to the UWP field, the algorithm now transfers 
to the solution of the global form of Eqs. 18 and 19. Use of a modified approach 
to Newton’s method is required during initial iterations (25). Entries in the 
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Jacobian A” are appropriately formed with v, or the ratio €/k held constant. 
The modified scheme and underrelaxation (a < 1) of the solution during initial 
iterations are required to prevent k and « from becoming negative. Computational 
experience shows that negative values of k and € quickly destroy the iteration 
process and cause the solution to diverge. 

A few modified iterations are required to create k and € values that are 
partially compatible with the current UWP values. The iterations are not continued 
until the k and ¢ solutions fully converge, because the current UWP solution 
is only approximate. Control is now passed back to the mean flow equation 
solver with updated values of k and « and therefore new values of v,. Several 
iterations are performed with a = | to align the UWP values with the updated 
distribution of v,. The solution scheme continues by passing from one equation 
set to another, iterating within each set as required, until there is no significant 
change in any of the variables. After several cycles the full Newton method 
may be implemented for the k-e equations. Element Jacobian terms for this 
solution process are listed in Appendix I. 


Bounpary Conpitions 


Figure 1 depicts a basin in which the fluid enters at the upper left corner 
and exits over a weir at the upper right corner. The top of the flow is bounded 


FIG. 2.—Finite Element Discretization of Tank Domain 


by a surface of constant (atmospheric) pressure called a free surface. The sides 
and bottom of the basin are bounded by solid impermeable walls. The bottom 
of the tank is usually inclined at some small angle 8. For a circular tank, 
the effect of slowly moving sludge rakes is ignored. The rakes in a rectangular 
basin often operate intermittently, in which case the rake velocity is periodically 
zero. 

A typical finite element discretization of a tank domain is shown in Fig. 
2. Along the inlet plane AB the values of U, W, k, and «€ are specified as 
essential boundary conditions. The value of W is generally specified to be zero. 
Values of U, k, and € should be obtained or inferred from experimental 
measurements. If such data are not available, these quantities may be estimated 
from other suitable computations or experiments. For example, the flow charac- 
teristics of a plane jet may aid the analyst in selecting appropriate inlet boundary 
values for this particular example. 

The exit region EFG is approximated in the numerical computations. Although 
the free surface in this region actually curves downward under the influence 


A F 
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of gravity as the fluid moves from F to G, the effect is of only local significance. 
A proper model of this region would substantially increase the overall computa- 
tional cost since many additional elements would be required. For this reason 
the exit boundary of the flow field is conveniently chosen to be FE, which 
is far enough upstream of the weir that the surface is nearly level. The exit 
flow is modeled as a point sink centered at G and of such strength as to match 
the inflow discharge. This approximation appears reasonable since the head 
on the weir in most tanks is never more than 1% or 2% of the tank depth. 
The velocity components along EF are therefore computed in accordance with 
a sink located at point G. As experimental data become available, actual measured 
velocities in this exit region could be used for the velocity boundary conditions. 
Boundary values of k and « for the exit region are undoubtedly determined 
from upstream events and the proximity of the free surface AF and vertical 
boundary DE. In this case it seems best to let the program compute the exit 
values of k and «. Normal derivatives of k and « along EF are thus set to 
zero. This specification is not strictly correct, since the flow from EF to G 
is nonuniform; however, any error introduced in this region will most likely 
not propagate far upstream because of the high velocities encountered in this 
region. If and when experimental values of k become available, essential boundary 
conditions along EF may then be specified, and the dissipation rate may be 
estimated from « = k°’?// in which | = a typical eddy length scale of the 
flow in this region. 

The free surface A F is modeled as a rigid lid (18). The computational boundary 
AF is fixed as a horizontal surface of zero slope. Pressures along this surface 
are computed however, and thus the effects of any small surface slopes are 
still accounted for, since the pressure along the free surface is allowed to vary 
[the term (dP/dr) in the horizontal momentum equation is nonzero]. The free 
surface boundary is also treated as a symmetry plane so that W = 0, the shear 
stress uw is zero so that (@ U/az) and (@ W/ax) = 0, and the normal derivatives 
(ak/az) and (d€/4z) are zero. 

The boundary conditions for the near-wall region are prescribed in a manner 
similar to those proposed by Launder and Spalding (17). The turbulence model 
is only valid in a fully turbulent regime. Therefore the finite element mesh 
is placed a distance 5 away from the solid boundaries (see Fig. 2) and beyond 
both the viscous sublayer and the buffer region. The velocity component normal 
to the wall is required to vanish. Along the vertical walls BC and DE this 
condition requires U = 0, while along CD, W = U tan B. A “‘slip’’ boundary 
condition is used parallel to the walls; the tangential stress is specified as a 


natural boundary condition. The wall shear stress is computed from the logarithmic 
law of the wall 


in which V = velocity parallel to the wall; x = von Karman constant (~0.4); 
U, = friction velocity; E = constant (~9 for a smooth wall); and y* = yU,R,. 
Here y = normal distance from the wall to the point in question and R, = 
(U,h,)/v. The value of y = 8 is selected to lie in the inertial subrange where 
the flow is assumed to be fully turbulent (y* > 30) but sufficiently close to 
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the wall (y ~ < 400) so that the shear stress is nearly constant (5). The shear-stress 


boundary condition along [, in Eqs. 16 and 17 is specified using Eq. 23 as 
follows: 


tar = -| r™N,c,U|U|1, dT 
r, 


r, 


| uwl,dT = — | r" N,c,W|W\1,dT 
re 
in which c,= x’ [In(Ey* )] The coefficient c ‘resembles a resistance coefficient 
and is a function of U,. 

In the near-wall region production nearly balances dissipation or A = e. 
Combined with Eqs. 7 and 23, this relation provides the following wall boundary 
conditions on k and e: 


k= 
€ = 


A pressure datum is required to obtain a unique solution. This is accomplished 
by setting P = constant at any one interior node and deleting the corresponding 
continuity equation from the global equation set for that node. The pressure 
may not be specified along any inlet or exit plane node, since conservation 
of mass along these boundaries then cannot be enforced (6). 


ApPLICATION—RECTANGULAR SEDIMENTATION BasiN 


The two-dimensional planar model is used here to determine the velocity 
field within a rectangular settling unit. The predictions are obtained for a basin 
of modest size; length = 40 ft (12.2 m), width = 15 ft (4.6 m), inlet depth 
= 10 ft (3.0 m), outlet depth = 9 ft (2.7 m). Recommended operational 
surface-loading rates (discharge entering basin / surface area of tank) for untreated 
wastewater range from 600-1,200 gal/day/sq ft (0.028-0.057 cm/s) (19, pp. 
448) and may be as high as 1,600 gal/day/sq ft (0.075 cm/s), (1, pp. 323) 
for plants with primary settling of iron- and polymer-coagulated raw wastewater. 
A surface loading rate of 1,500 gal/day/sq ft (0.071 cm/s) was selected for 
this example. This rate should simulate conditions in an overloaded primary 
clarifier and may indeed represent the lower limit of clarifier performance. 
For this example the characteristic velocity and length scale are chosen to be 
U, = 0.1 ft/sec (0.0305 m/s) and h, = 10 ft (3.05 m). 

The finite element grid used in the numerical simulation is presented in Fig. 
2. The distribution of elements is selected to be dense in the inlet and exit 
regions (where velocity gradients are large) and relatively more sparse in the 
interior domain. 

The inlet region for this problem was chosen to be the upper three elements 
on the left hand side of the computational grid; the inflow velocity profile 
was approximated as one half of a plane jet. As noted by Larsen (14) in his 
full-scale studies of rectangular settling units, simple inlet structures produce 
velocity patterns similar to free jets a short distance from the inlet. A literature 
search for values of U, k, and € pertaining to a plane jet (2,10) aided in selecting 
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approximate inlet conditions. The remaining boundary conditions were described 
in the previous action. 

A vector plot of the computed velocity field is depicted in Fig. 3. The jet 
expands to fill the tank approximately 2.5 tank depths downstream from the 
inlet. The presence of the outlet is felt thereafter as the fluid streamlines converge 
toward the upper right hand portion of the basin. A recirculation zone is formed 
below the expanding jet. The zeros plotted in the lower corners represent velocities 
which are too small to be plotted by the graphic routine. In fact, close inspection 
of the computer output revealed a small recirculating eddy in each of these 
corners. At this time, experimental data of detailed point velocity measurements 
in primary sedimentation basins are unavailable for comparison. It is encouraging 
to note, however, that flow patterns measured by Clements and Khattab (3) 


= 

=> 
= = 

= 
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FIG. 4.—Contour Plot of Dimensionless Turbulent Kinematic Viscosity v, 


on scale-model circular basins with surface inlets are qualitatively similar to 
the computed flow pattern depicted in Figure 3. 

A contour plot of the dimensionless turbulent kinematic eddy viscosity [v, 
= v*/(h, U,)] is shown in Fig. 4. The distribution of v, is clearly not uniform; 
the most pronounced changes occur in the inlet and exit regions where velocity 
gradients are largest. Near the solid boundaries v, must decrease due to the 
damping effect of the wall on the turbulent fluctuations. Elsewhere the variation 
is less pronounced. It is quite apparent that this distribution of v, would be 
very difficult to specify a priori. 

The numbers of elements and node points for this example are 135 and 454 
respectively. The numbers of mean flow equations (UWP) is 995, and the number 
of turbulence equations (ke) is 788. The approximate numbers of nonzero entries 
in the Jacobian matrix A is 114,000 and 76,000, respectively, for the UWP 
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and ke sets of equations. This example required 6.2 min of central processor 
time on a CDC 7600 computer. The total number of UWP and ke iterations 
required to obtain a fully converged solution is 53. 

A useful application of this solution is to use the known velocity field to 
solve the now-linear convection-diffusion equation describing the transport of 
particulate matter (25). In fact, the predicted variation of v, also allows one 
to assign known nonconstant values to the transport coefficients appearing in 
this equation. The settling efficiency of the basin may then easily be determined 
from the computed distribution of particulate concentration along the boundary 
of the tank. 


Summary ano Conciusions 


A mathematical model for turbulent flow in primary sedimentation basins 
has been formulated and solved numerically via the Galerkin finite element 
method. Because the eddy viscosity is determined as part of the overall solution, 
the uncertainty about the selection of values for these coefficients is resolved. 
The predictive abilities of the model are therefore enhanced. 

The computed velocity field clearly indicates the existence of a nonuniform 
recirculating flow pattern within the basin. These results qualitatively confirm 
the experimental measurements of several studies (3,14,28). The computed 
distribution of v, varies considerably throughout the settling basin. 

This project has shown that the modeling of flow in settling basins can be 
improved in several ways by additional research. First, a series of detailed 
experimental studies which can supply point velocity and possible Reynolds 
stress measurements in full-scale settling basins are needed to allow one to 
assess properly the validity of the current model. Once validated, the computa- 
tional algorithm could be used to improve settling performance by examining 
the effects of tank geometry or inlet and exit configurations on predicted removal 
rates, or both. Second, while the present method can successfully compute 
the solution to a complex recirculating turbulent flow, the sensitivity of the 
Newton-based solution algorithm to initial values for the variables seems 
excessive. Additional research should endeavor to develop a solution procedure 
whose performance is less dependent on the initialization. Finally, the current 
model should be augmented to include stratified phenomena. Either temperature 
or species concentration variation may be the causative mechanism. However, 
mathematical models of turbulent stratified flow phenomena are less well 
developed than those for simpler flows and will undoubtedly require improvement 
and allied experimental confirmation before such flows are well simulated. 


ACKNOWLEDGMENT 


This paper is based on work supported by the National Science Foundation 
under Grant No. ENG 7618846. 


Aprenoix 1.—Evement Entries in Jacosian Matrix 


aU aN aN 
wl u—+ wins) + da 
au, a, or ar az 


SEDIMENTATION BASINS 


N, aN, k 
0z 


ON, aN, 
— +—-—— dn 
or 
0k aN, ak an, 
or ar 


HY5 587 
(.8N, aN, aN, an, 
+ —_-—+] an 
a, or oar 
aN, auw 
or au, 
4 au aN, aN 
r™N,N,—d2+\ r™v,—— ao 
aw, a, az a, ar 
aN, 
or 
4 an 
mN,M,dQ—\ r™M,—dO 
oP, a, a, or 
4 aw aN, aN 
r™N,N,—dQ+\ r"v,——da 
au, a, or a, or 
aN, 
r, 
aN ow oN oN, aN 
aw, a, or a, or ar 
aN, aN, auw an, 
+2——J dn -|\ Jar ....... 
az az r aw, az 
aN 
r™M,—dQ2+\ G3) 
oP, a, az r, 
a a 
au, a, or 
aw, a, 
4 4 
+ 
a, 
c 
a, 


MAY 1981 


k? 
| +1] dQ 
a, € 


+ 


2 
€ 

-| + aa 
a, k 


de ON, de ON, 

4+ 
or or az 

aN, aN, 2e 
+ 

or oz 
at aN,aN, aN, aN, 


o.\ or ar 


-| ar" N,N,c,|U\1, aT 


auw 
=-\ 2"N,N,c,|W\1, aT 
r, aw, 


r, 


Appenoix 


. Wastewater Treatment Plant Design, Manual No. 36, by a Joint Committee of the 
American Society of Civil Engineers and the Water Pollution Control Federation, 
ASCE, 1977. 

. Bradbury, L. J. S., “‘The Structure of a Self-Preserving Turbulent Plane Jet,”’ Journal 
of Fluid Mechanics, London, England, Vol. 23, 1965, pp. 31-64. 

. Clements, M. S., and Khattab, A. F. M., ‘‘Research into Time Ratio in Radial Flow 
Sedimentation Tanks,’’ Proceedings of the Institution of Civil Engineers, Vol. 40, 
Aug., 1968, pp. 471-494. 

. Gibson, M. M., and Launder, B. E., ‘“‘On the Calculation of Horizontal, Turbulent, 
Free Shear Flows Under Gravitational Influence,”’ Journal of Heat Transfer, American 
Society of Mechanical Engineers, Feb., 1976, pp. 81-87. 

. Gosman, A. D., and Ideriah, F. J. K., “‘Teach-T: A General Computer Program 
for Two-Dimensional, Turbulent, Recirculating Flows,’’ Report of the Department 
of Mechanical Engineering, Imperial College, London, England, 1976. 

. Gresho, P. M., Lee, R. L., Sani, R. L., and Stullich, T. W., ‘‘On the Time-Dependent 
FEM Solution of the Incompressible Navier-Stokes Equations in Two- and Three- 
Dimensions,’’ UCRL-81323, Lawrence Livermore Laboratory, July, 1978. 

. Gresho, P. M., and Lee, R. L., ‘Don’t Suppress the Wiggles—They’re Telling You 
Something!"’ Proceedings of the Winter Annual American Society of Mechanical 

' Engineers Meeting, Dec., 1979. 

. Ha Minh, H., and Chassaing, P., ‘“‘Some Numerical Predictions of Incompressible 
Turbulent Flow,’’ Numerical Methods in Laminar and Turbulent Flows, C. Taylor, 
et al., eds., Pentech Press, London, England, 1978, pp. 287-300. 


588 HY5 
_ 
a, 
af 
ak, 
| 
af. 
af. ) an 
de ON, ON, 
a, oO, € or ar 
a uw 
aU, 
1 
2 
3 
4 
5 
6 


HY5 SEDIMENTATION BASINS 589 


9. 


10. 


Hamlin, M. J., and Wahab, A. H. A., ‘‘Settling Characteristics of Sewage in Density 
Currents,”’ Water Research, Vol. 4, 1970, pp. 609-626. 

Hanjalic, K., and Launder, B. E., ‘“‘A Reynolds Stress Model of Turbulence and 
its Application to Thin Shear Flows,” Journal of Fluid Mechanics, London, England, 
Vol. 52, 1972, pp. 609-638. 


. Hanjalic, K., and Launder, B. E., ‘“‘Contribution towards a Reynolds-Stress Closure 


for Low-Reynolds-Number Turbulence,’’ Journal of Fluid Mechanics, London, Eng- 
land, Vol. 74, 1976, pp. 593-610. 


. Ideriah, F. J. K., “On Turbulent Forced Convection in a Square Cavity,’’ Numerical 


Methods in Laminar and Turbulent Flow, C. Taylor, et al., eds., Pentech Press, 
London, England, 1978, pp. 257-269. 


. King, I. P., Norton, W. R., and Iceman, K. R., “‘A Finite Element Model for 


Two-Dimensional Flow,’’ Finite Element Methods in Flow Problems, J. T. Oden, 
et al., eds., UAH Press, 1974, pp. 133-137. 


. Larsen, P., “On the Hydraulics of Rectangular Settling Basins Experimental and 


Theoretical Studies,’’ Report No. 1001, Department of Water Resources Engineering 
Lund Institute of Technology, University of Lund, Lund, Sweden, 1977. 


. Launder, B. E., Morse, A., Rodi, W., and Spalding, D. B., ‘‘Prediction of Free 


Shear Flows—A Comparison of the Performance of Six Turbulence Models,”’ Proceed- 
ings of the Conference on Free Turbulent Shear Flows, NASA Langley Research 
Center, Vol. 1, July, 1972, pp. 361-422. 


. Launder, B. E., Reece, G. J., and Rodi, W., ‘‘Progress in the Development of a 


Reynolds-Stress Turbulence Closure,’’ Journal of Fluid Mechanics, London, England, 
Vol. 68, 1975, pp. 537-566. 


. Launder, B. E., and Spalding, D. B., ‘“‘The Numerical Computation of Turbulent 


Flows,’’ Computer Methods in Applied Mechanics and Engineering, Vol. 3, 1974, 
pp. 269-289. 


. McGuirk, J., and Rodi, W., ‘“‘A Depth-Averaged Mathematical Model for the Near 


Field of Side Discharges into Open-Channel Flow,”’ Journal of Fluid Mechanics, 
London, England, Vol. 86, 1978, pp. 761-781. 


. Wastewater Engineering: Collection Treatment, Disposal, Metcalf and Eddy, Inc., 


McGraw-Hill Book Co., Inc., New York, N.Y., 1972. 


. Montens, A., ‘“‘The Use of Radioactive Isotopes for Water Flow and Velocity 


Measurement,”’ Radioisotope Conference 1954, J. E. Johnston, ed., Butterworths, 
London, England, pp. 169-180. 


. Oliver, A. J., “‘A Finite Difference Solution for Turbulent Flow and Heat Transfer 


Over a Backward Facing Step in an Annular Duct,’’ Numerical Methods in Laminar 
and Turbulent Flow, C. Taylor, et al., eds., Pentech Press, London, England, 1978, 
pp. 467-478. 


. Patankar, S. V., and Spalding, D. B., Heat and Mass Transfer in Boundary Layers, 


2nd ed., Intertext, 1970. 


. Sani, R. L., Gresho, P. M., and Lee, R. L., ‘On the Spurious Pressures Generated 


by Certain GFEM Solutions of the Incompressible Navier-Stokes Equations,’’ Third 
International Conference on Finite Elements in Flow Problems, June, 1980. 


. Schamber, D. R., and Larock, B. E., ‘“‘A Finite Element Model of Turbulent Flow 


in Primary Sedimentation Basins,’’ Proceedings of the Second International Conference 
on Finite Elements in Water Resources, C. A. Brebbia, et al., eds., Pentech Press, 
London, England, 1978, pp. 3.3-3.21. 


. Schamber, D. R., ‘‘Finite Element Analysis of Flow in Sedimentation Basins,’’ thesis 


presented to the University of California at Davis, in 1979, in partial fulfillment of 
the requirements for the degree of Doctor of Philosophy. 


. Schamber, D. R., and Larock, B. E., ““Computational Aspects of Modeling Turbulent 


Flows by Finite Elements,’’ Computer Method in Fluids, C. Taylor, et al., eds., 
Pentech Press, London, England, 1980, pp. 339-361. 


. Schamber, D. R., and Larock, B. E., ‘Constant and Variable Turbulent Eddy Viscosity 


Flow Simulations in Rectangular Settling Tanks,’’ Presented at the Aug. 6-8, 1980, 
ASCE Hydraulics Division Specialty Conference, held at Chicago, Ill. 


. Wills, R. F., and Davis, C., ‘‘Flow Patterns in a Rectangular Sewage Sedimentation 


Tank,”’ Presented at the 1962, International Conference on Water Pollution Research, 
held at London, England, pp. 335-385. 


_ 


590 MAY 1981 
Aprenoix Ill.—Notation 


The following symbols are used in this paper: 


Jacobian of equation matrix; 

resistance coefficient; 

constants in turbulence model; 

turbulent diffusion of quantity in parenthesis; 
constant in law of wall relation; 

Froude number; 

equation residuals; 

system of nonlinear equations; 

ratio of weight to mass; 

elevation of free surface; 

reference length; 

turbulence kinetic energy; 

typical eddy length scale; 

direction cosines; 

basis functions for interpolating polynomials; 
constant defining planar or axisymmetric flow; 
iteration counter; 

mean pressure; 

production term; 

(U,h,)/v (Reynolds number); 

radial / horizontal coordinate; 

mean radial /horizontal velocity component; 
reference velocity; 

friction velocity; 

radial / horizontal turbulent velocity fluctuation; 
velocity parallel to solid boundary; 

rake velocity; 

circumferential turbulent velocity fluctuation; 
mean vertical velocity component; 

vertical turbulent velocity fluctuation; 
horizontal coordinate; 

solution vector to f; 

corrections applied to x; 

dimensionless normal distance from wall boundary of 
tank; 

yU,Ro; 

relaxation parameter; 

angle of inclination of basin floor with horizontal; 
fluid domain boundary; 

element boundary; 

distance from basin wall to computational grid; 
rate of dissipation of turbulence kinetic energy; 
von Karman constant; 

molecular kinematic viscosity; 

turbulent kinematic viscosity; 
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fluid density; 
pUuUu, pvv, pww, puw turbulent Reynolds stresses; 
o,,0, constants in turbulence model; 
® velocity gradient contribution to production term; 
2 fluid domain; and 
2, element domain. 
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UNCERTAINTIES RESULTING FROM CHANGES 
IN RIVER ForRM 


By Durl E. Burkham' 


INTRODUCTION 


Two recent laws enacted by the United States Congress have focused attention 
on the need for more information on the processes of sedimentation and on 
the hydraulics of flood flow in alluvial channels. The two laws are the Federal 
Water Pollution Control Act Amendments of 1972 (Public Law 92-500) (10) 
and the Flood Disaster Protection Act of 1973 (Public Law 93-234) (11). Public 
Law 92-500 requires, among other things, that an environmental impact statement 
for the affected area be made before actions are taken to manage the products 
of a drainage basin—water and debris—and to manage flood plains in the basin. 
Public Law 93-234 requires, among other things, that flood-prone areas be 
delineated. These requirements, which place man in the context of nature, demand 
that the reasons for changes in the landscape and the side effects of the changes 
be fully understood. 

Major channel-form changes in a drainage basin can affect all components—the 
channel, surface-water yields and movement, sediment yields and movement, 
ground-water supply and movement—of the drainage system and the interrelation 
of the components. As will be analyzed in more detail subsequently, a specific 
river-form change in a large drainage basin may originally be local; however 
the effects of a local change usually propagate along alluvium-filled valleys 
to other parts of the basin. 

This paper is concerned primarily with hydrologic implications and uncertainties 
that are created by temporal changes in alluvial-channel form in relatively large 
watersheds. Historical evidence of river-form changes is presented first. In 
sequential sections, brief examinations are presented about hydrologic implica- 
tions and uncertainties relative to channel-form changes and land use, flood 
characteristics of alluvial streams, and hydrology and sediment yields for basins 
in the Southwest. Changes in river form are produced by processes involving 
long-term geologic influences, annual and seasonal fluctuations in the water 
and sediment runoff, and artificial influences exerted on the river by man’s 


activities. This paper is particularly concerned with the aforementioned short-term 
processes. 


"Hydro., United States Geological Survey, 2800 Cottage Way, Sacramento, Calif. 95825. 

Note.—Discussion open until October 1, 1981. To extend the closing date one month, 
a written request must be filed with the Manager of Technical and Professional Publications, 
ASCE. Manuscript was submitted for review for possible publication on December 19, 
1978. This paper is part of the Journal of the Hydraulics Division, Proceedings of the 
American Society of Civil Engineers, ©ASCE, Vol. 107, No. HY5, May, 1981. ISSN 
0044-796X /81 /0005-0593 / $01.00. 
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The term ‘“‘uncertain,”’ as used in this paper, means not definite, ascertainable, 
or fixed, as in time of occurrence, number, quantity, or place. Uncertainty— 
unpredictability, indeterminacy, or indefiniteness—is manifested as errors in 
data and in predictions. 

River-channel form (river form) is the plan-view pattern and vertical cross-sec- 
tional shape of the river channel between river banks. Alluvial river form is 
described qualitatively as meandering, straight, or braided, and narrow and deep, 
or wide and shallow. Changes in river form cause changes in the flood plain. 
The flood plain is defined as nearly level land that occupies the bottom of 
the valley of a present stream and is subject to flooding unless protected artificially 
(1). A flood, as used in this report, is defined as the occurrence of water 
in excess of channel capacity; overbank flow occurs which inundates parts 
or all of the flood plain (15). For this paper, anything that occurs naturally 
on a flood plain—e.g., alluvial fans at the mouths of tributaries or vegetation—is 
considered part of the flood plain. 

River-channel form, although a small part of the total landscape, is of special 
interest to man because development traditionally has been along the fertile, 
level land of the flood plains of rivers. Changes in the river-channel form, 
therefore, affect the quality of life and the value of property for these human 
intruders on the flood plain. 

Short-term processes effect changes in river form in alluvial channels in all 
regions of the United States. For most rivers, especially those in nonmountainous 
humid climates, the changes are usually minor, and the hydrologic consequences 
of the changes are insignificant. The changes in many rivers, especially some 
of those in the semiarid Southwest, however, have created highly significant 
consequences. Evidence of significant river-form changes for a few rivers in 
different locations in the United States are given in this paper. The factors 
and mechanisms involved in causing these changes are examined only in a 
general way. 


Historica Evioence of River-Form Cuance 


The Rio Salado in central New Mexico is one of many streams in the Southwest 
where changes have occurred. The Rio Salado, a tributary of the Rio Grande, 
rises on the north side of the Datil Mountains and has a general eastward 
course north of that range to a junction with the Rio Grande at the village 
of San Acacia. According to Bryan (7), when Lorenzo Padilla, an early settler, 
saw the valley along the Rio Salado in 1880, the channel near Santa Rita, N.M., 
was not large and the broad, flat valley seemed a propitious place for farming. 
Daniel Curry, while surveying near Santa Rita in 1882, recorded the width of 
the stream bed as ranging from about 12 ft-49 ft (4 m-15 m) along a number 
of section lines (7). 

Another surveyor in 1918, according to Bryan (7): 


. . . found the course of the river radically different from that shown 
in Curry’s survey of 1882, his measurements ranging from 330 to 550 
ft in the same stretch of stream channel where Curry found widths of 
11.88 to 48.84 ft. 


HY5 RIVER-FORM CHANGES 
Bryan (7) further states that: 


. . . unlike many similar streams in New Mexico, which have not only 
widened their channels but deepened them in the same period, the Rio 
Salado, at least in the vicinity of Santa Rita, has even yet banks that 
are only 3 to 10 ft high and average about 5 ft high. It is obvious, however, 
that the whole regimen of the stream is much different from that which 
existed in 1880. 


Smith (18) reported that the Republican River, in the northwestern corner 
of Kansas, was, prior to 1935, a narrow stream with a practically perennial 
flow of clear water, and with well-wooded banks. The Republican River in 
1940, however, had a broad, shallow, sandy channel with intermittent flow; 
the trees had been destroyed, much valuable farmland on the valley bottom 
was sanded over, and the channel had been filled several feet. Smith associated 
the river-form change with the occurrence of a major flood in 1935. 

Schumm and Lichty (16) reconstructed the river history of the Cimarron River 
in southwestern Kansas for the period 1874-1960. They reported that, in 1874, 
the river averaged 50 ft (15 m) in width in six counties. In the period 1914-1942, 
however, the Cimarron River widened to about 1,200 ft (366 m). After 1942, 
reconstruction of the Flood plain was in progress. By 1954, the channel averaged 
about 550 ft (168 m) in width. 

Hack and Goodlett (12) reported that a violent cloudburst flood of June 1949 
caused severe erosion in the Little River Valley in the mountain region in the 
Central Appalachians of Virginia. They state: 


The runoff produced dozens of debris avalanches on the upper mountain 
slopes, enlarged most of the channel ways, and reworked the debris on 
the botton lands of many larger valleys and in places removed the forest 
cover on the entire valley floor. 


A rare flood in December 1964 caused river-form changes in several valleys 
in northern California. Stewart and LaMarche (19) wrote: 


The December 1964 flood on Coffee Creek (a small high-gradient mountain 
stream in Trinity County, California) was of rare frequency and unprece- 
dented in historic time. Erosion and deposition during the flood were 
catastrophic and significantly changed the character of the valley. Some 
of this erosion and deposition was similar to that described by Hack and 
Goodlett (12) for northern Virginia. Within the valley, the preflood channel 
was commonly filled, and new channels formed at entirely different 
locations. 


Storms in 1969 caused dramatic changes along small valleys in urban southern 
California. According to Scott (17): 


A unique combination of substantial channel change and documentation 
of the changes by high-order photogrammetry was studied in Tujunga 
Wash in southern California. Extensive scour and fill occurred during 
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the record breaking 1969 floods in this 3-mile long, partly urbanized fanhead 
valley. Maxima of 20 + 2 ft of net scour and 35 + 2 ft of net fill were 
measured. . . . Net elevation change of the channel thalweg varied from 
as much as 14 + 2 ft of scour to as much as 16 + 2 ft of fill. 


The writer (4) studied changes in river form in the Gila River in Safford 
Valley in southeastern Arizona for the period 1846-1970. The changes can be 
grouped into three periods—1846-1904, 1905-1917, and 1918-1970. From 1846- 
1904, the stream channel was narrow and meandered through a flood plain 
covered with willow, cottonwood, and mesquite. Only moderate changes occurred 
in the width and sinuosity of the stream channel in this period; the average 
width of the stream channel was less than 150 ft (46 m) in 1875 and less than 
300 ft (91 m) in 1903. 

The average width of the stream channel of the Gila River increased during 
1905-1917 to about 2,000 ft (610 m), mainly as a result of large winter floods 
that carried relatively small sediment loads (4). The meander pattern of the 
stream and the vegetation in the flood plain were destroyed completely by 
the floods. 

The stream channel of the Gila River narrowed during 1918-1970 and the 
average width was less than 200 ft (61 m) in 1964 (4). The stream channel 
developed a meander pattern, and the flood plain became densely covered with 
vegetation. Minor widening of the stream channel occurred in 1965 and in 1967, 
and the average width of the channel was about 400 ft (122 m) in 1968. The 
reconstruction of the flood plain in 1918-1970 was accomplished almost entirely 
by the accretion of sediment. 


The most important factors influencing the deposition of sediment in the 
Gila River during 1918-1970 were the wide flood channel and the small floods 
that carried relatively large sediment loads (4). The large sediment loads resulted 
mainly from rapid erosion of the alluvial deposits in the drainage basins tributary 
to the Gila River. 


Hyorotoeic Impuications AND Uncertainties Due to Cuances in River Form 


Introduction.—Major floods apparently were a primary cause of the river-form 
change for each of the samples cited in the preceding section. Large flood 
flows exert great force on the channel banks and on objects in the main flow 
path. During a major flood, the main flow path often is straight down valley. 
While the meander pattern is intact, part of the flow is directed along the 
meandering stream channel and a large amount of turbulence develops along 
the streambanks. As a result of the stresses produced by the turbulent forces 
along the streambanks and around other stationary objects, channel changes 
eventually take place—banks erode, trees are uprooted and flushed downstream, 
grass is removed, and alluvial fans at the mouths of tributaries are destroyed. 
These changes result in a wider and cleaner stream channel that is more conducive 
to rapid movement of flood waves and transport of the basin’s debris during 
floods. 

Reconstruction of a flood plain, typically during periods of relatively low 
flood peaks, is accomplished naturally by accretion of sediment; accretion results 
because the river, in absence of major floods, does not have the hydraulic 
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competence to move the inflow of sediment farther downstream. Flood plains 
may be built by sediment accretion in five general ways: (1) By the development 
of islands in the stream channel and their subsequent attachment to one bank 
by channel abandonment; (2) by direct deposition on the flood plains; (3) by 
deposition in the stream channel along the banks; (4) by formation of natural 
levees; and (5) by deposition on alluvial fans at the mouths of tributary streams. 
Rapid reconstruction of a flood plain depends on rapid erosion of surfaces 
or channels in the upstream drainage basin. 

To generalize, much is known about the natural erosion and sedimentation 
processes involved when river-form changes occur. Most textbooks dealing with 
sediment transport and geomorphology present adequate general descriptions 
of the factors and processes involved. Based on present knowledge, many general 
questions concerned with the movement of water and sediment can be resolved. 
Public Laws 92-500 and 93-234, however, encompass specific problems pertinent 
to the movement of water and erosional debris that often require precise answers 
to specific questions. Often, exact determinations of the effects of man’s activity 
are sought. Precise answers to specific questions, however, are not always 
obtainable, especially when the results of recent channel-form changes caused 
by nature are involved. For example, according to Wolman (20) 


. . Tuling on a controversy over logging in a redwood forest in California, 
Judge R. H. Kroninger wrote the following: ‘‘While numerous expert 
witnesses in the field of geology, forestry, engineering, and biology were 
presented, their conclusions and the opinions they derived from them 
are hopelessly irreconcilable in such critical questions as how much and 
how far solid particles will be moved by any given flow of surface water. 
They were able to agree only that sediment will not be transported upstream”’ 
[State of California, Marin County versus E. Richetti, and others, 1969]. 


As will be examined subsequently, the controversy over logging involves channel 
changes caused by floods. 

Related reasons why flood-induced channel changes cause uncertainties perti- 
nent to the aims of Public Laws 92-500 (10) and 93-234 (11) are as follows: 


1. The capability to adequately determine whether a river is susceptible to 
a river-form change has not been established. 

2. The capability to predict when a river-form change will occur does not 
exist because of the preceding reason and because of an inability to predict 
when a channel-changing major flood will occur. 

3. A specific river-form change in a large drainage basin may originally be 
local; however the effects of a local change usually propagate along alluvium-filled 
valleys to other parts of the basin (4,9). The propagation process is not understood 
except in a general way (20); the rate of propagation through specific parts 
of a basin is varied and the rate cannot be adequately predicted. 

4. The capability to predict the rate of reconstruction of a flood plain does 
not exist because, among other things, the rate depends on the sequence of 
floods, on flow rates, and on the sequence and rates of sediment inflow. 

5. An adequate accounting of the timing and discontinuities of processes 
involved in sedimentation presently is not possible because the physical laws 
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of sediment transport are incompletely known, especially for the complex problem 
of transportation of the wide range of sediment sizes that are carried by natural 
streams, and sedimentation in natural streams involves many interrelated vari- 
ables, most of which cannot be assigned meaningful numerical values. 

6. For the first four reasons, drainage basin changes resulting from man’s 
activities cannot be determined adequately because a basis for identifying 
flood-induced modification of a basin—river-form changes and subsequent 
adjustments—has not been established. 


Because of the six aforementioned reasons, many problems pertinent to the 
management of water and debris, the management of flood plains, and the 
effects of man’s activities in a drainage basin cannot be resolved except in 
a probabilistic or general way. The analysis in subsequent sections relates 
uncertainties caused by flood-induced channel changes to specific hydrologic 
problems and to specific procedures used in hydrologic studies. 

Land-Use Problems.—Public Laws 92-500 (10) and 93-234 (11), as previously 
implied, require that the effects of man’s use of drainage basins be fully 
understood. The controversies and uncertainties involving logging in the redwood 
forest in California and overgrazing in the Southwest are examples of those 
that involve uses of basins where flood-induced channel changes have occurred. 

Logging in Redwood Forest in California.—The controversy in the redwood 
forest in California involves the effects of logging in causing changes in the 
forested landscape of the downslope Redwood National Park. The harvesting 
of the native timber probably influences changes in water discharge and sediment 
load, as well as adjustments in width, depth, sinuosity, gradient, and location 
of stream channels in the forest. Concurrently, natural flood-flow modification 
unrelated to the timber harvesting, if any, is likely to have been, or will be, 


of the same type, although not necessarily in the same direction. Janda (14) 
states: 


Considerable geologic and geobotanical evidence is accruing to indicate 
that devastating floods (that is, large floods that bring about major change 
in stream channel and hill slope configuration) like those that occurred 
in 1955, 1964, and 1972 are naturally occurring phenomena in this area. 
. . . The impact of these floods on stream sediment load is both immediate 
and of long duration... . 


A long-term result of the large floods is to increase the amount of readily 
erodible sediment in and immediately adjacent to the stream channel. 

Another probable long-term consequence of large floods in the redwood forest 
that must be accounted for in resolving the controversies over logging involves 
the translation of the effects of river-form changes to other parts of the drainage 
system. As stated by Gilbert in Chorley, Dunn, and Beckinsale (9): 


Of the main conditions which determine the rate of erosion, namely; quantity 
of running water, vegetation, texture of rock, and declivity, only the last 
is reciprocally determined by rate of erosion. . . . Wherever by reason 
of change in any of the conditions the erosive agents come to have locally 
exceptional power, that power is steadily diminished by the reaction of 
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rate of erosion upon declivity. Every slope is a member of a series, receiving 
the water and the waste of the slope above it, and discharging its own 
water and waste upon the slope below. If one member of the series is 
eroded with exceptional rapidity, two things immediately result: first, the 
member above has its level of discharge lowered, and its rate of erosion 
is thereby increased; and second, the member below, being clogged by 
an exceptional load of detritus, has its rate of erosion diminished. The 
acceleration above and the retardation below, diminish the declivity of 
the member in which the disturbance originated; and as the declivity is 
reduced the rate of erosion is likewise reduced. 


Precise answers such as those sought by the presiding judge (in the controversy 


over logging) presently are not obtainable because of the reasons previously 
given. Wolman (20) states: 


Judge R. H. Kroninger’s pique was probably directed at the fact that 
he had to choose a single number as the proper criteria for legal action 


(in the controversy over logging) and no good scientist would or could 
give him a single number. 


Overgrazing in Southwest.—The controversy concerning the influence of 
extensive grazing on erosion in the Southwest in the latter part of the nineteenth 
century centers around the relative effects of land use as compared to natural 
flood-flow modification. Conservationists, for many years, have promoted the 
theory that overgrazing was a primary cause of the rapid erosion that occurred 
in alluvial valleys in the Southwest during the period 1880-1930. The period 
when rapid erosion gutted valleys followed closely years when large numbers 
of livestock were brought into the area and the amounts of precipitation were 
small. Logically, the combination of little precipitation and extensive grazing 
should cause a deterioration in the vegetation of the valleys, which would have 
made the alluvium more susceptible to erosion. Considering only these factors, 
the theory that overgrazing was a primary cause of deterioration of alluvium-filled 
valleys is easily accepted. Considerable evidence, however, is now available 
to indicate that overgrazing, in fact, may have been only a minor factor. 

The beginning of severe erosion in southwestern valleys apparently correlates 
exceptionally well with the incidence of channel-changing major floods. Major 
floods, therefore, may have been a primary factor in causing the severe erosion 
in most of the southwestern valleys. 

Bryan (7) stated: 


Valid conclusions as to the merits of these theories [theories concerned 
with the causes of severe erosion in the Southwest] cannot be reached 
until historical data on the time at which erosion began have been 
accumulated. Knowledge of the date of the beginning and progress of 
this spectacular change in the regimen of streams is particularly necessary 
in arriving at a decision as to the effects on erosion processes of the 
introduction of cattle and sheep and the overgrazing that in most localities 
ensued. 
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He associated the beginning of severe erosion along the Rio Salado with an 
exceptional rain and flood in 1883. The Rio Salado drains a basin where extensive 
grazing occurred in the latter part of the nineteenth century. Bryan (7) stated: 


. . . the progress of erosion [in the Rio Salado] begun near Santa Rita 
in 1883 has been, as measured in years, fairly slow and has not yet, 
43 years later, affected all the minor tributaries. 


The writer’s observations during a visit to the site in 1969 indicate that head-cut 
erosion is still progressing upstream in minor tributaries; this probably indicates 
that the cycle of erosion that began with the flood in 1883 is continuing. In 
1969, the depth and width of the Rio Salado near Santa Rita apparently were 
not significantly different than those recorded in 1918. 

The results of a study by Schumm and Lichty (16) seem to indicate that 
overgrazing was an insignificant factor, and floods may have been a primary 
factor in causing the channel changes along the Cimarron River in southwestern 
Kansas. They wrote: 


The period of channel widening was initiated by the maximum flood of 
record, which was followed by a long period of deficient precipitation, 
and was terminated by the major flood of 1942. The period of flood-plain 
construction was characterized by above-average annual precipitation and 
floods of low to moderate peak discharge.. . . McLaughlin. . . concluded 
that, because there was no major change in number of livestock, grazing 


was not an important factor determining channel changes. . . . Livestock 
density was greatest during the period of channel narrowing 1942-48. The 
activities of man and animals in the Cimarron River basin undoubtedly 
had their effects on the regime of the river, but any relation between 
cultivation, grazing, and channel change is obscure. 


Studies of channel changes along the Gila River in southeastern Arizona 
apparently also support a premise that floods were the major factor in causing 
channel changes in the Southwest and that overgrazing probably was a minor 
factor (4). The widening of the Gila River in 1891, 1905-1917, 1941, and 1965-1967 
was coincident with major floods. Most of the floods originated in the mountainous 
part of the headwater area as a result of large storms. The floods originated 
in areas where there was no significant grazing. Grazing apparently did not 
have a significant influence in causing the major floods and probably had only 
a minor influence in causing the widening of the stream channel of the Gila 
River. Large-scale grazing began in about 1872 in parts of the Gila River below 
the shaded mountain forest and below the area that produced most of the flood 
flow; grass was the dominant vegetation in these areas. If grass had been the 
dominant vegetation in the flood-producing areas, then the floods probably would 
have been slightly more severe. 

The high flows in the Gila River drainage during the period 1905-1917 lowered 
the local base level and accelerated erosion in the tributaries of the Gila River 
by increasing the gradient of most of the tributaries at their confluence with 
the Gila River. High flows in the tributary streams provided the motive power 
necessary to start erosion (4). Generally, severe erosion occurred rapidly in 
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the steep streams that drain the mountainous terrain near the Gila River; however 
the areal extent of the channel erosion was small because of the small areal 
extent of the easily erodible alluvium that underlies these streams. Overgrazing 
was not a factor in causing erosion in the mountainous terrain because grazing 
was insignificant. Overgrazing could have been a contributing factor in causing 
erosion in alluvial valleys drained by gently sloping streams because they were 
grazed extensively, and erosion in some of the valleys was severe. 

The San Simon River, tributary to the Gila River in Safford Valley, is an 
example of a gently sloping stream that has undergone severe channel erosion 
since 1905 (4). The San Simon River drains an area of about 2,200 sq mile 
(5,698 km”), and its valley covers most of the drainage basin. The San Simon 
River was an insignificant and poorly-defined watercourse in 1903. Debris from 
side tributaries apparently had been collecting in the alluvium-filled valley for 
centuries, and in 1905 it was poorly drained and relatively unstable. When severe 
erosion began, probably during the major floods of 1905, deep channels were 
cut and eventually became large; erosion then spread to the side tributaries 
according to Gilbert’s theory (9). Around 1919, a deep eroded channel had 
extended upstream for about 60 mile (97 km). By 1960, there were more than 
100 mile (161 km) of gullied channel from 10 ft-40 ft (3 m-12 m) deep, and 
from 20 ft-500 ft (6 m-152 m) wide. The gullied channels captured runoff that, 
prior to the erosion, would have spread over the valley and replenished soil 
moisture necessary for plant growth. As the water flowed into the deep channels, 
additional erosion occurred; however the eroded sediment, generally of small 
particle size, was easily moved downstream because the flow was confined. 
A ditch dug at the mouth of the San Simon River may have influenced channel 
erosion in the stream. The ditch was dug by settlers prior to 1900 to divert 
floods in the San Simon River away from the cultivated land. Because severe 
erosion took place at the same time in other gently sloping streams, the erosion 
in the San Simon River basin probably would have occurred even if the ditch 
had not been dug. Because severe erosion occurred in steep channels that drain 
the mountainous terrain near the Gila River where grazing was not a factor, 
the erosion in the San Simon River basin probably would have occurred even 
if extensive grazing had not taken place. 

To summarize, major floods that caused changes in channel form probably 
were a primary factor influencing severe erosion in alluvial valleys in the 
Southwest during 1880-1930. Extensive grazing during the latter part of the 
nineteenth century, however, probably was a contributing factor. Even after 
an elapsed time of about 100 yr and the collection of data during many studies, 
precise determinations of the effects of extensive grazing on erosion in relative 
large drainage basins in the Southwest have not been made primarily because 
of the reasons previously described. 

Flood Characteristics of Alluvial Streams 

Introduction.—Major channel-form changes can significantly affect the flood 
characteristics of streams—timing, transformation of floods in travel, magnitude 
and frequency of floods, and the bankfull capacity of streams. As will be indicated 
in the following, studies—such as those of flood routing, frequency of flooding 
and channel capacity, and in-regime flow—that deal with flood characteristics 
of streams, must consider the effects of channel-form changes in order to minimize 
errors and bias. 
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Flood Routing.—Flood routing is the process of determining progressively 
the timing and magnitude of a flood wave at successive points along a river. 
The prediction of the timing and magnitude of a flood wave in a reach of 
interest is usually based on a flood-routing model. The principal types of 
flood-wave movement are classified as uniformly progressive and reservoir action 
(8). Uniformly-progressive flood-wave movement is the downstream movement 
of a flood wave that does not change in shape; such movement occurs under 
ideal conditions in a prismatic channel in which the variation in resistance along 
the channel is small. Reservoir action refers to the transformation of a flood 
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of Flood Waves Moving Through a 55-mile (88-km) Reach of Gila River in Safford 
Valley, Ariz., for Indicated Water Years [after (5)) 


wave caused by reservoir pondage which results when the channel is irregular 
and the variation in resistance along the channel is large. Both types of flood-wave 
movement can occur at different times in the same reach; the type that will 
occur apparently is dependent largely on the river form. 

More than 250 inflow and outflow hydrographs for flood waves in 1914-1970 
for a 55-mile (88-km) reach along the Gila River in southeastern Arizona were 
reviewed to determine temporal adjustments in flood-wave transformation caused 
by changes in the river form (5). As previously noted, the channel was relatively 


40; 4 
30} 
1914. } 
40: 
30- 4 
1914-27 
20: 
| 
| 20 
40 4 
+ 
20 
~~ 
1914-97 if 
%/55~~-- 


HY5 RIVER-FORM CHANGES 603 


wide—about 2,000 ft (610 m)—in 1905-1917. During 1918-1940; the channel 
changed; the stream channel became narrow, stream channel meander increased, 
natural levees developed along the stream channel, flood plain vegetation spread 
and became dense, and large alluvial fans developed at the mouths of tributaries 
(4). After the changes, the Gila River consisted of a stream channel in which 
flow below bankfull discharge—3,000 cu ft/sec-5,000 cu ft/sec (85 m*/s-142 
m’/s) after 1930—moved relatively fast, and a congested flood plain in which 
flood flow moved slowly. Floods having average peak discharges of 10,000 
cu ft/sec-20,000 cu ft/sec (283 m*/s-566 m*/s) (mean of peak rates at the 
ends of the reach) moved through the 55-mile (88-km) reach in about 14 h 
in 1914-1927; however it took them more than 40 h to move through the reach 
in 1961-1970 (see Fig. 1). 

Flood waves during 1914-1927 generally retained their inflow shape (hydrograph 
shape at the upstream end of the study reach) as they moved through the 
study reach on the Gila River (see Fig. 2), and the flood-wave movement 
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FIG. 2.—Measured and Synthesized Flood Flow, July 14-16, 1919, and September 
3-5, 1925, at End of 55-mile (88-km) Reach of Gila River in Safford Valley, Ariz. 
[after (5)] 


approximated a uniformly progressive one (5). The uniformly-progressive flood- 
wave movement, however, applies only to floods that were not reduced signifi- 
cantly by infiltration. Infiltration during many floods in the 1914-1927 period 
significantly reduced the size of the floods, and in many instances, the inflow 
shapes were altered. 

Floods were transformed greatly from about 1935-1970 as they moved through 
the study reach on the Gila River, probably as a result of reservoir action 
caused by the meandering stream channel and flood plain vegetation (see Figs. 
3 and 4) (5). Infiltration also may have been a cause for the change in inflow 
shape. Significant temporal changes in the attenuation of peak rates moving 
through the study reach accompanied changes in wave shape. Peak rates of 
less than about 13,000 cu ft/sec (368 m’/s) were reduced to bankfull discharge, 
between 3,000 cu ft/sec and 5,000 cu ft/sec (85 m*/s and 142 m’/s), for 
most floods that occurred in 1955-1965, and the average reduction in peak 
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rates ranged from about 7,000 cu ft/sec (198 m’/s) for floods having inflow 
peaks of 14,000 cu ft/sec (396 m’*/s) to about 4,000 cu ft/sec (113 m’/s) 
for floods having inflow peaks of 44,000 cu ft/sec (1,250 m*/s) during the 
same period. The increased attenuation of peak flow from 1914-1927 to 1944-1965 
probably resulted from temporal increases in reservoir action and infiltration. 
Decreases in the amount of streamflow contributed by tributary watersheds 
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FIG. 3.—Measured and Synthesized Flood Flow, July 23-25, 1955, at End of a 55-mile 
(88-km) Reach of Gila River in Safford Valley, Ariz. [after (5)] 


along the study reach may also have caused some difference in peak outflow. 
Tributary streamflow ponded behind natural levees during 1944-1965 (5); the 
natural levees were not present during 1914-1927. 

Conclusions resulting from the study on the 55-mile (88-km) reach on the 
Gila River that may be pertinent to flood routing in other rivers subject to 
river-form changes are: 


1. The size and meander pattern of the stream channel of the Gila River 
are determined by past dominant flows. The stream channel is wide and straight 
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at the end of a period in which high flows were dominant, and is narrow and 
has a meander pattern at the end of a period in which low flows were dominant. 

2. The stream channel and flood plain system, when fully developed for a 
dominant flow, has a persistent effect on floods. A low-flow system—developed 
by and for low flows—attenuates flood peaks passing through the reach; the 
peak flows of flashy floods may be reduced to bankfull discharge. A high-flow 
system—developed by and for high flows—does not increase flood rates; however 
streamflow from side tributaries along the study reach may contribute more 
significantly to peak rates in the Gila River when a high-flow rather than a 
low-flow system exists. 

3. The downstream velocity of the center of mass of flood waves that had 
peak discharges of between 10,000 cu ft/sec and 20,000 cu ft/sec (283 m*/s-566 
m’/s) during 1914-1927 may have been as much as three times that for the 
same rates during 1943-1970. 
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FIG. 4.—Measured and Synthesized Flood Flow, January 11-19, 1960, at End of 
a 55-mile (88-km) Reach of Gila River in Safford Valley, Ariz. [After (5)] 


4. A low-flow system may change rapidly to a high-flow system when a 
series of major floods occurs; however several years of low flow are required 
before a high-flow system changes to a low-flow system. It took about 50 yr 
for the 1970 low-flow system in the Gila River to develop. 

5. Annual peak flows measured at the downstream end of the study reach 
reflect, among other things, the persistent effect of the upstream system, and 
therefore they are not as random in time as are flows at the upstream end. 
Because of changes in the system, the data of peak flows collected at the 
downstream end of the study reach during 1914-1927 are from a different 
population than the data of peak flows for the period 1943-1970. 

6. A flood-routing model used to predict the timing and transformation of 
a flood wave must have the capability to account for the effects of river-form 
changes; this accounting is difficult to accomplish for many streams because 
of reasons previously given. 
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Frequency of Flooding and Channel Capacity.—The frequency of flooding 
in a reach of river clearly depends on both the flow rate and the bankfull 
capacity of the stream channel. River-form changes along a river, therefore, 
can affect the frequency of flooding in two ways—by changing the peak flow 
rate and by changing the bankfull capacity of the river. This causes uncertainties 
in the delineation of flood-prone areas. 

The method of delineating flood-prone areas—flood mapping—for flood-in- 
surance studies pertinent to Public Law 93-234 usually consists of two parts: 
(1) Peak discharge-frequency analysis; and (2) hydraulic analysis. A peak 
discharge-frequency curve is a product of the frequency analysis. This curve 
is based mainly on historical peak discharge rates; however precipitation data 
are used in conjunction with annual peak flows for some analyses. These rates 
reflect, among other things, the persistent effects of the upstream drainage 
system and the effects of river-form changes in the system. As previously 
mentioned, the data for a given site for a flood occurring during a period when 
a stream channel is narrow and meandering through a valley are often from 
a different population than data for a flood occurring in a period when the 
channel is wide and straight. As shown by the writer (5), the reduction in 
peak flow caused by a major channel change can be significant. 

Hydraulic analyses needed to predict 7-year water-surface elevations (the 
elevation that will occur, on an average, once in T-years—10 yr, 50 yr, 100 
yt) for flood-insurance purposes generally assume that the present channel 
boundary will not change significantly prior to the occurrence of the 7-year 
flood. This assumption is not valid for many streams embedded in unconsolidated 
alluvial deposits, and it may lead to significant bias in predicting T-year elevations. 
The standard error of estimate for T-year elevations is relatively large even 
for rigid boundary conditions; the nationwide standard error of estimate is about 
5.0 ft (0.9 m) (6). If the channel boundary is not rigid, then major changes 
in river form can result when extreme floods occur, and the accuracy for the 
predicted water-surface elevation probably is low even though the magnitude 
of discharge for a given frequency of occurrence is assumed to be predictable. 
For example, river-form changes along the Gila River that have been previously 
described probably resulted in a 6 ft-8 ft (1.8 m-2.4 m) difference in the water 
surface for a 25-yr flood. 

In-Regime Concepts.—Many engineering approaches to solving hydraulic 
problems in alluvial channels use the in-regime concepts introduced by Blench 
(2). Blench’s definition of ‘‘in-regime”’ is: 


. . that average values of the quantities we appreciate as constituting 
regime do not show a definite trend over some interval—usually of the 
order of a score or two of years . . . [rivers in regime] demonstrate 
themselves to us in the form of varying discharges, breadths, depths, 
velocities, meander patterns, sediment contents, and so forth, but their 


average behavior does not usually change greatly over small periods of 
historic time. 


The in-regime approach to solving hydraulic problems has been to describe 
the hydraulic characteristics of a river in terms of time-average quantities. The 
hydraulic geometry equations developed by Leopold and Maddock (15) are 
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examples of this approach. The in-regime concept is valid for streams where 
the channel boundary can withstand stresses to such an extent channel form 
is not significantly changed during temporary periods of high flow. Even though 
the boundary of a stream embedded in unconsolidated alluvium may be in a 
constant state of change when flow occurs, the in-regime concept may be valid 
if major channel-form changes do not occur. 

The major channel-form changes for alluvial valleys described previously in 
this paper have occurred mainly because of temporal changes in the dominant 
flow rate. Small (relative to a channel that could contain a major flood) channels 
that have a meander pattern and a vegetation-covered flood plain are closely 
associated with dominantly low-flow rates. The in-regime concept probably applies 
to these relatively small channels and low-flow conditions (relative to rates 
during a major flood). The in-regime concept also probably applies to relatively 
large straight channel and high-flow conditions. Considering a long period of 
time and a large region where channel-form changes occur, most channels are 
small because of the temporal and spatial dominance of low-flow rates. For 
the period and region, the in-regime concept probably is valid for most channels. 
A few channels in the large region, however, may be in a state of flux between 
a high-flow channel and a low-flow channel; the in-regime concept is not valid 
for this condition and these channels. Of the channels in a state of flux, some 
stay relatively large for a long period of time even though low-flow conditions 
prevail. The Rio Salado in New Mexico is an example. Of the channels in 
a state of flux, some have changed in the last 100 yr from a low-flow type 
to a high-flow type and back to a low-flow type; the Gila River in Safford 
Valley is an example. 

In-regime concepts are used as a basis for estimating average runoff from 
drainage basins in arid and semi-arid regions (13). In the procedure, measured 
average annual flow is regressed against selected dimensions of channel ge- 
ometry—usually width and depth at bankfull stage; or width and depth of a 
cross section between bars and berms. The degree of correlation usually is 
reasonably good because, on an average, most channels reflect the temporal 
dominance of low-flow rates in a region and the average annual discharge 
adequately represents the low-flow rate. The regression equation is used to 
estimate runoff from basins drained by ungaged streams. The procedure usually 
gives adequate results if data for channels that are in a state of flux between 
a high-flow channel and a low-flow channel do not significantly affect results. 

Hydrology and Sediment Yields for Basins in Southwest.—Major river-form 
changes can cause significant differences in the surface runoff, in recharge 
to aquifers, and in sediment yields in affected basins, especially in arid and 
semi-arid regions of the United States. Most of the natural recharge to aquifers 
in arid and semi-arid regions occurs as a result of infiltration during intermittent 
flow in channels embedded in alluvium deposits. Infiltration depletes runoff 
and replenishes ground water. Infiltration in an alluvial channel is dependent 
on many factors (3): (1) Surface area and depth of water in the region; (2) 
permeability, moisture distribution, and temperature of the subsurface alluvium; 
(3) physical quality of the water and length of time the water is available at 
the land surface; (4) chemical quality of the surface and subsurface water; 
and (5) structural stability of the porous media. Many of the dominant forces 
and parameters affecting infiltration are related directly to the width, depth, 
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and velocity of streamflow. As previously shown, these parameters are signifi- 
cantly altered by major changes in channel form. A significant change in the 
infiltration rate would cause a significant change in the runoff and the ground-water 
recharge. 

The rate of sediment movement at a site in a stream depends on two 
factors—hydraulic conditions at the site and the supply of sediment readily 
available for movement. Major channel-form changes significantly alter both 
factors. The water-discharge to sediment-discharge relationship for a site, 
therefore, often is significantly transformed when a major channel-form change 
occurs; the annual sediment yield from the affected basin also usually is 
significantly altered. 


Summary ano ConcLusions 


The uncertainties resulting from major changes in river form in many large 
drainage basins are of such magnitude that the objectives of Public Laws 92-500 
and 93-234 relative to the management of water and debris, the management 
of flood plains, and the effects of man’s activities presently cannot be met 
except in a probabilistic or general way. Rivers in arid and semi-arid regions 
are more susceptible to major changes in river form than those in more humid 
regions. 

Channel-form changes caused by natural processes create uncertainties for 
the following related reasons: 


1. The capability to adequately determine whether a river is susceptible to 
a river-form change has not been established. 

2. The capability to predict when a river-form change will occur does not 
exist because of the preceding reason, and because of an inability to predict 
when a channel-changing major flood will occur. 

3. Aspecific river-form change ina large basin may originally be local; however 
the effects of a local change usually propagate along alluvium-filled valleys 
to other parts of the basin. The propagation process is not understood except 
in a general way, and the rate of propagation through specific parts of a basin 
is varied and cannot be adequately predicted. 

4. The capability to predict the rate of reconstruction of a flood plain does 
not exist. Among other things, rate depends on the sequence of floods, on 
flow rates, and on the sequence and rates of sediment inflow. 

5. An adequate accounting of the timing and discontinuities of processes that 
are involved in sedimentation presently is not possible because the physical 
laws of sediment transport are incompletely known, especially for the complex 
problem of transporting the wide range of sediment sizes that are carried by 
natural streams, and sedimentation in natural streams involves many interrelated 
variables, most of which cannot be assigned meaningful values. 

6. For the first four reasons, water basin changes resulting from man’s activities 
cannot be determined adequately because a basis for identifying changes in 


a drainage system resulting from natural flood-flow modifications has not been 
established. 


Examples of controversies and perplexing situations that involve channel-form 
changes are: 
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1. Logging in the uplands of the redwood forest in California: The controversy 
involves the role of logging in causing changes in the forested landscape of 
the downslope Redwood National Park. The harvesting of the native timber 
probably influences changes in the rates of erosion and deposition in the forested 
land and influences changes in width, depth, sinuosity, gradient, and location 
of stream channels in the forest. Concurrently, natural flood-flow modification 
unrelated to timber harvesting is likely to have been, or will be, of the same 
type, although not necessarily in the same direction. Because of reasons outlined 
in the preceding, changes in the basin caused by man’s activity cannot be precisely 
established. 

2. Overgrazing in the Southwest: Conservationists, for many years, have 
promoted the theory that overgrazing was a primary cause of the rapid erosion 
in alluvium-filled valleys that occurred in the Southwest during the period 
1880-1930. Considerable evidence is now available to indicate that overgrazing, 
in fact, may have been only a minor factor. The rapid erosion apparently was 
a naturally-occurring phenomenon brought about largely by major floods. 

3. Flood routing: Prediction of the timing and transformation of a flood wave 
based on a flood-routing model can have significant error if an accounting of 
the effects of river-form changes has not been properly considered. 

4. Frequency of flooding and hydraulics of floods: The frequency of flooding 
in a reach of river depends on both the flow rate and the bankfull capacity 
of the stream channel. River-form changes affect the frequency of flooding 
in two ways—by changing the peak flow rates and by changing the bankfull 
capacity of the stream channel. Flood-insurance studies pertinent to Public Law 
93-234 are made using the assumption that the present bounding conditions 
existing in a channel will not change significantly prior to the occurrence of 
a T-year flood. This assumption can lead to significant error and bias for channels 
where major channel-form changes have occurred. 

5. In-regime concepts: The in-regime concept is valid for streams where the 
channel boundary can withstand stresses to such an extent that channel form 
is not significantly changed during high flows. Even though the boundary of 
a stream embedded in unconsolidated alluvium may be in a constant state of 
change during a variation of flow rates, the in-regime concept is valid if major 
channel-form changes do not occur. Many channel-form changes result because 
of temporal changes in the dominant flow rate. Relatively small channels correlate 
with relatively low flow rates; relatively large channels correlate with relatively 
high flow rates. Many channels in a large region may be in a state of flux 
between a high-flow channel and a low-flow channel; the in-regime concept 
is not valid for these channels. 

6. Surface runoff, ground-water recharge, and sediment yields: Major river 
changes can cause significant differences in the surface runoff, recharge to 
aquifers, and sediment yields in the affected basins. 
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MAXIMUM CLEAR-WATER SCOUR 
AROUND CIRCULAR PIERS 


By Subhash C. Jain," M. ASCE 


INTRODUCTION 


The safe and economical design of bridge piers requires accurate prediction 
of the maximum expected depths of scour of the stream bed around them. 
The interaction between the flow around a bridge pier and the erodible sediment 
bed surrounding it is very complex. In fact, the phenomenon is so involved 
that only very limited success has been enjoyed by the attempts to model scour 
computationally, and physical models remain the principal tool employed for 
estimating expected depths of scour. Two types of scour may be identified: 
(1) Clear-water scour—where material is removed from the scour hole and not 
replaced; and (2) scour—that occurs with general sediment transport. Following 
the experimental study of Chabert and Engeldinger (8) on local scour around 
bridge piers, most of the investigators in this area concurred on the general 
shape of the curve which delineates the variation of scour depth with mean 
flow velocity. According to this curve, scour depth increases with increase 
in mean velocity in the clear-water regime, reaches an absolute maximum at 
a velocity approximately equal to the threshold velocity (hereinafter referred 
to as the mean velocity for incipient sediment motion), and decreases slightly 
with further increase in mean velocity in the sediment-transport regime where 
it fluctuates nonperiodically about the equilibrium scour depth due to bed-form 
migration. Since the maximum scour depth is required in designing bridge piers, 
most of the experimental studies in the past were conducted either in the 
clear-water regime or with flow velocities not much higher than the threshold 
velocity in the sediment-transport regime. 

A wide variety of empirical equations based upon a limited range of data 
(both laboratory and prototype) have been developed in the past to estimate 
the maximum scour depths around bridge piers. Unfortunately the relatively 
large scatter in the available data on local scour around bridge piers makes 
it possible to fit a wide variety of curves which diverge greatly at high Froude 
numbers and high relative depths of the flow. Not only that their extrapolation 
to higher Froude numbers and relative depths cannot be used as a sound basis 
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for pier design, it is also difficult to draw conclusions regarding the appropriateness 
of these formulas at low Froude numbers and relative depths. In this study, 
the potential predictors of the maximum clear-water scour are compared with 
the experimental data. The limitations of these predictors are considered. Another 
formula to predict the maximum clear-water scour is proposed. 


Estimation oF Scour Dertu 


The flow in the vicinity of the pier is so complex that a complete analytical 
or numerical description of the scour process is not possible at the present 
time. Accordingly, phenomena involving local scour around bridge piers have 


TABLE 1.—Constants and Exponents in Eq. 1 


Breusers (6) 


incipient 


rational 


motion 

Larras (16) incipient rational 
motion 

Blench (5) sediment regime* 
transport 

Laursen sediment rational 
transport 

Laursen and Toch (19)°| incipient rational 
motion 

Arunachalam (3) sediment regime" 
transport 

Ahmad (1) sediment regime 
transport 

Shen, et al. (28) clear water rational 

Shen, et al. (28) incipient rational 
motion 

Hancu (12) incipient rational 
motion 

Inglis-Poona incipient rational 


(Thomas, (29)] motion 


“y = regime depth. 
*Transformed and simplified by Melville (20). 


“Transformed and simplified by Neill (27). 
*K = 1.2. 


been studied most extensively in laboratory experiments, from which several 
empirical formulas have been developed to estimate the maximum scour depths 
around bridge piers. In general, they are based upon a limited range of data 
and are applicable to conditions similar to those for which they were derived. 
It is difficult to confirm their adequacy for design purposes due to limited 
field measurements. Though there are some similarities among the various 
empirical relations, they differ widely in terms of the hydraulic variables 
considered to be significant. Most of the scour relations can be expressed in 
the form of the following general equation 


(1) (2) (3) (4) (5) (7) (10) 
I 1.4267") 0} 0 | O | O 
Il 1.80 -1}3/4] 0 | 1] 0 
Il Ll 
1.35 0}03|/0 0 
ul 1.95 -1|5/6} 0 | 1] 0 
Ill 3.18K* |-1| 1 [2/3] 1 | 0 
Ill 11.00 C1 rt 
3.40 0 }1/3|2/3| 0 | 
2.42 0 |1/3|2/3| 0 | 
4.05 1] 0 


CLEAR-WATER SCOUR 


in which d, = the scour depth measured below mean bed elevation; y = the 
flow depth; F = V/V gy = the Froude number; V = the flow velocity; 5 
= the width of the pier projected on a plane normal to undisturbed flow; D 


TABLE 2.—Scour Relations not Expressible in Form of Eq. 1 


Ap- 
Group] Investigator | Regime | proach Remarks 
(3) (4) (6) 
Breusers, et d,/b = 2 tanh (y/b) 
al. (7) 
Chitale (9) i d,/y =(-0.51 
+ 6.65F — 5.49F7) 
Inglis Lacey D, = 0.946 (Q/f)'”” D,, = scour depth below 
(13) water surface 
Q = discharge, in cubic 
feet per second 
f = Lacey silt factor 
= 1.76 V 
Knezvic (15) rational | d, = 8.72 [(q q=y 
Bata (4) rational | d,/y = 10 [F? 
3D/y) 
Maza (20) rational | d,/b =f (F,y/b) graphical form 


Hancu (12) rational | d,/b V.. = threshold velocity 


= 2.42(2V/V, gy 
— (y/d)'” 
Garde (11) D,/y a =(B-—b)/B 
= 4.0 4, 7273 1/a(F)" | B= clear channel width 
Ns and n 
are functions of particle 
drag coefficient, Froude 
number, and pier shape 
Chabert and |clear rational | d, =f (b,y, V,D) graphical form 
Engel- water 
dinger (8) and sed- 
iment 
transport 


= the mean sediment size; A and B = constants; m, n, p, and r = exponents. 
The values of the constants A and B, and the exponents m, n, p, and r depend 
upon the pier shape, the angle of attack of the flow, and the sediment properties. 
Eq. | indicates that the scour depth is a function of four variables; the pier 
size, the flow depth, the flow velocity, and the sediment size. The values of 
the constants and exponents in Eq. | for circular piers based on the relations 
proposed by the various investigators are summarized in Table 1. There are 
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other empirical relations which could not be expressed in the form of Eq. 1; 
these are listed in Table 2. 

The scour relations can be classified in several ways. They are grouped into 
six categories depending upon the number of significant hydraulic parameters 
in each relation. The only significant parameter is the pier size in group I; 
the flow depth and the pier size in group II; the flow depth and velocity in 
group III; the flow depth and velocity, and the sediment size in group IV; 
the flow depth and velocity, and the pier size in group V; all four parameters 
in group VI. The scour formulas are additionally classified into three groups 
based on the flow regimes; (1) Clear water; (2) incipient sediment motion; and 
(3) sediment transport regimes. The formulas for the condition of the incipient 
sediment motion predict the maximum clear-water scour which is assumed to 
occur at flow velocity approximately equal to the threshold velocity; therefore 
F in Eq. | should be replaced by the threshold Froude number, F.(F.=V./V gy, 
in which V, = the threshold velocity). The value of F in Eq. 1 for clear-water 
scour should be less than F.. The empirical relations are categorized into two 
classes based on two different approaches; (1) Regime; and (2) rational approaches. 

A comparison made by Anderson (2) of the several formulas listed in Tables 
1 and 2 showed that the estimates of the relative scour depths differed widely, 
particularly for the higher values of Froude number and relative depth. The 
divergence among the curves representing the various scour formulas clearly 
indicated that most of these equations are applicable only for certain range 
of flow conditions and should not be extrapolated for flow conditions outside 
that range. A summary of the experimental data including the range of parameters 
used to develop the various scour formulas is presented elsewhere (14). 


Maximum Crear-Water Scour 


The dimensional analysis of the scouring parameters has been presented by 
several investigators including Breusers, et al. (7) and Neill (23). On the assumption 
that the influence of fluid viscosity on scour is negligible, the scour depth 
for natural sediments can be expressed as 


Furthermore, the relative sediment size, D/y, can be expressed as a function 
of threshold Froude number F. using the logarithmic velocity distribution and 
the Shields’ criterion for the initiation of sediment movement (see the curve 
in Fig. 2.46 of Ref. 30, labeled Shields), i.e.: 


D 


Eq. 2, after replacing D/y by F. from Eq. 3, reduces to 


d, y 


The maximum clear-water scour occurs at F = F.; Eq. 4 for F = F. can be 
simplified to 


CLEAR-WATER SCOUR 


TABLE 3.—Summary of Data Used in Comparative Analysis 


Flow Mean Maximum 
: velocity, Pier sediment Flow scour Thresh- 
in meters |diameter,| size, in depth, depth, old 
per in centi- milli- in centi- | in centi- | Froude 
second meters | meters meters | meters | number® Investigator 


(1) (7) 


The writer and 
Fischer 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Shen, et al. 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
Chabert and 
Engeldinger 
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0.82 5.1 2.50 24.7 8.7 0.46 
0.40 15.2 0.24 21.9 18.0 0.21 
0.32 15.2 0.24 11.6 13.4 0.27 
0.36 15.2 0.24 15.6 15.8 0.24 
0.38 15.2 0.24 20.6 17.1 0.22 
0.44 15.2 0.24 21.0 21.0 0.22 
0.41 15.2 0.24 26.3 18.6 0.20 
0.38 15.2 0.46 17.6 16.6 0.25 
0.50 91.4 0.46 61.0 54.9 0.15 
0.85 5.0 3.00 20.0 7.0 0.56 
0.85 10.0 3.00 20.0 12.5 0.56 
0.85 15.0 3.00 20.0 18.5 0.56 
0.76 5.0 3.00 10.0 8.7 0.71 
0.76 10.0 3.00 10.0 13.1 0.71 . 
0.76 15.0 3.00 10.0 17.5 0.71 . 
0.66 5.0 1.50 20.0 9.8 0.39 
0.66 10.0 1.50 20.0 17.0 0.39 
0.66 15.0 1.50 20.0 20.3 0.39 
0.40 5.0 0.52 19.7 9.5 0.25 
- 0.40 10.0 0.52 19.7 12.2 0.25 
0.40 15.0 0.52 19.7 14.9 0.25 
0.42 5.0 0.52 35.0 9.0 0.20 
0.42 10.0 0.52 35.0 12.0 0.20 
0.42 15.0 0.52 35.0 13.7 0.20 
0.37 10.0 0.52 10.0 11.5 0.32 ; 
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TABLE 3.—Continued 


Chabert and 
Engeldinger 
Hancu 
“The value of F. are based on the threshold velocity from the Shield’s criterion for 
the critical shear stress and the logarithmic velocity distribution given by V./V, = 2.5 
In (11.02 yx/D,,), when V, = the critical shear velocity, and x = a correction factor 


accounting the effect of viscosity (30). A value of 0.01 cm?/s for the kinematic viscosity 
of water was assumed in the computation. 


Shen et al (28) 


Chabert & Engeidinger (8) 


Join & Fischer (14) 


Hancu (12) 


Breusers (Correlation coefficient,r= ?) 
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FIG. 1.—Comparison of Breusers’ Formula with Available Scour Data 
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CLEAR-WATER SCOUR 


Only those predictors among the various formulas listed in Tables 1 and 2, 
which are in the form of Eq. 5, and based on rational approach, were selected 
for the comparative analysis. 

Scour formulas based on “‘regime’’ approach used an empirical formula, derived 
mainly from canal data, relating the average stable channel depth of an alluvial 
channel to the dominant discharge and bed material. Some scour formula assumed 
that scour depths at structures could be expressed as some multiple of the 
average regime depth. These equations should be applied in cases in which 
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Loursen (r= 0.75) 
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FIG. 2.—Comparison of Laursen’s Formula with Available Data 


the flow, sediment transport, and channel characteristics are quite similar to 
those from which the particular formula was derived. The equations based on 
“‘regime’’ approach are, therefore, not considered in the comparative analysis. 
Furthermore, the scour formulas which predict scour depth independent of pier 
size are applicable only for small values of the relative depth. These formulas, 
such as in groups III and IV in Tables | and 2, also are excluded from the 
comparative analysis due to their limited applications. Chabert and Engeldinger 
(8) did not propose any predictor for scour depth. Their experimental data 
are utilized in the analysis. The scour formulas of group V in Tables | and 
2 are in the form of Eq. 5. This equation would be in accord with the scour 
formulas of group II if the influence of F. on scour depth is found to be small. 
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FIG. 3.—Comparison of Hancu’s Formula with Available Scour Data 
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FIG. 4.—Comparison of Formula by Shen, et al. with Available Scour Data 
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If the effect of the depth of flow in addition is found to be insignificant, Eq. 
5 would belong to the scour formula of group I. In order to determine the 
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FIG. 6.—Comparison of Maza’s Formula with Available Scour Data 
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FIG. 7.—Comparison of Garde’s Formula with Available Scour Data 


range of flow parameters for which the various scour formulas are valid, these 
formulas are compared with some of the available scour data for circular piers. 
The experimental data of Chabert and Engeldinger (8), Shen, et al. (28), Hancu 
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(12), and the writer and Fischer (14) were used in the comparative analysis 
because all tests in these studies were conducted with circular cylinders and 
all basic experimental data were easily available. Only flow conditions with 
0.02 < (F — F.) < 0.1 were included. A summary of experimental data in 
the comparative study is presented in Table 3. The lower limit of 0.02 instead 
of zero was used to insure that the equilibrium scour depth was achieved in 
the experiment in a reasonable long period; the scour depth in clear-water regime 
approaches a limit asymptotically and takes a long time to reach this limit. 

The comparison between the observed scour depths and the scour depth 
predicted by various formulas is shown in Figs. 1 through 11. The values of 
linear correlation coefficient, r, between the observed and predicted scour depth 
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FIG. 8—Comparison of Formula by Larras with Available Scour Data 


are also included in these figures. The correlation coefficient indicates how 
well the equation fits the data. The lines of perfect agreement based on formulas 
by Larras, Breusers, et al., Shen, et al., and Laursen and Toch form an envelope 
for all data. Breusers’ relation (Fig. 1) envelops the data for (y/b) =< 3 only 
and thus underpredicts scour depths for ( y /b) = 3. The scour formula by Laursen 
is similar to that by Laursen and Toch as the former was derived by adjusting 
a constant in the solution for long constriction to agree with the latter. Laursen’s 
relation (Fig. 2), however, underpredicts for y/b < 1. The expressions by Shen, 
et al., and Hancu are identical except for the coefficient A (Eq. 1). Hancu’s 
formula-(Fig. 3) with A = 2.42 seems to give the best fit for the data while 
Shen’s formula (Fig. 4) with A = 3.40 envelops the data. Inglis-Poona’s relation 
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(Fig. 5) is in agreement with data for fine sand (D < 0.5 mm) and y/b =< 
4. Maza graphical relation agrees well with the data (Fig. 6) except it does 
not envelop the data of Shen, et al., and Hancu for fine sand and low values 
of y/b. The scour formula by Garde underpredicts for fine sand and low values 
of y/b and overpredicts for coarse sand and high values of y/b, as shown 
in Fig. 7. 

From design and safety considerations a predictor which envelops the experi- 
mental data is desirable. The formulas by Larras (Fig. 8), Breusers, et al. (Fig. 
9), Shen, et al. (Fig. 4), and Laursen and Toch (Fig. 10) fall in this category. 
A regression analysis of the experimental data (14), based on the assumption 
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FIG. 9.—Comparison of Formula by Breusers, et al. with Available Scour Data 


that the relative scour depth is a power function of the threshold Froude number 
and the relative depth, yielded 
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The coefficient 1.41 in Eq. 6 was increased by 30% to 1.84 in order to form 
an envelope for all data; it resulted in 
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The comparison of Eqs. 6 and 7 with the data is presented in Fig. 11. 

In order to determine the range of flow parameters for which the various 
formulas that envelope the experimental data, predict satisfactorily i.e., less 
overpredict, the following criterion was used. A formula was considered satisfac- 
tory if the difference between the predicted and observed values is less than 
30% of the predicted value. The value of 30% was chosen because the constant 
1.41 in Eq. 6 was increased by 30% to envelope the data. This criterion is 
represented by a dashed line in Figs. 4 and 8-11. This criterion is met by 
the formula of Breusers, et al. for coarse sand, Larras for y/b = 3, Shen, 
et al. for fine sand and y/b = 1, and Laursen and Toch for almost all of 
the data. It can, therefore, be inferred that the scour formula of Laursen and 
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FIG. 10.—Comparison of Formula by Laursen and Toch with Available Scour Data 


Toch and Eq. 7 are the best among these formulas to predict the maximum 
clear-water scour. The exponent of (y/b) in Eq. 7 is identical to that in the 
Laursen and Toch formula. These two formulas predict the same scour depth 
for F. = 0.29. Eq. 7 for F. < 0.29 predicts less scour depths than that given 
by Laursen and Toch formula. The following considerations show that Eq. 
-7 is a better predictor for the maximum clear-water scour than the Laursen 
and Toch formula. The laboratory experiments on local scour conducted recently 
by Ettema (10), Nicollet (24), and Hancu (12) indicate that the maximum 
clear-water scour depends upon the sediment size. Substituting for F, from 
Eq. 3 in Eq. 7 shows that scour depth is a function of sediment size, while 
Laursen and Toch equation predicts that it is independent of sediment size. 
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The correlation coefficient for Eq. 7 is higher than that for Laursen and Toch 
relation. For most practical cases F. is small and probably less than 0.3; Eq. 


7 predicts less (therefore leading to an economical design) but safe (because 
it forms the envelope for all data) scour depth. 


Conc.usions 


Several predictors for the maximum clear-water scour depth were compared 
with the available experimental data. The range of flow parameters for which 
these formulas either overpredicted or underpredicted were delineated. Predictors 
which enveloped all data were identified. The comparison indicated that the 
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FIG. 11.—Comparison of Eqs. 6 and 7 with Available Scour Data 


scour formula by Laursen and Toch (19) is the best predictor among those 
compared in this study, as it envelopes all data and overpredicts less than 
the other formulas. However, the Laursen and Toch formula predicts that the 
scour depth is independent of sediment size. Another formula (Eq. 7) for the 
maximum clear-water scour is proposed, which is very similar to that of Laursen 
and Toch but includes the effect of sediment size on scour depth. 
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Il.—Notation 


The following symbols are used in this paper: 


A,B constants; 
pier size; 
mean sediment diameter; 
scour depth; 
V/V gy = Froude number; 
” IV gy = threshold Froude number; 
gravitational constant; 
exponents; 
mean flow velocity; 
threshold velocity; and 
flow depth. 
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TECHNICAL NOTES 


Note.—Discussion open until October 1, 1981. To extend the closing date one month, 
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ASCE. This paper is part of the Journal of the Hydraulics Division, Proceedings of 
the American Society of Civil Engineers, OASCE, Vol. 107, No. HY5, May, 1981. 
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BIOFILM GROWTH AND HyDRAULIC PERFORMANCE 


By Brian Butterworth’ 


Rheidol is a 56-MW hydroelectric station operated by the Central Electricity 
Generating Board (of England) in mid Wales; its scheme is shown in Fig. 1. 
Both the low-pressure NANT-Y-MQCH to DINAS and the higher-pressure 
DINAS to CWM RHEIDOL tunnels have suffered biofouling and its attendant 
loss of hydraulic performance. 

The problem first manifested itself in the higher pressure tunnel, as a steady 
decline in the maximum power output. Twelve years after commissioning, this 
deficit had reached 10%, and was showing signs of a steadily increasing rate 
of deterioration. 

Detailed performance testing of the plant showed the turbines and generators 
to be performing as designed, while the pipeline exhibited a frictional head 
loss some 44% greater than design. 

Internal inspection of the pipeline revealed a layer of material some 4.5 mm 
thick, which was extremely hard and knobbly when wet (see Fig. 2); the material 
became soft and friable when dry. Analysis of the deposit showed a very high 
iron and manganese content, as well as organic material: SiO,—30.9%; loss 
on ignition at 800° C—30.2%; Zn—0.5%; Fe—6.9%; Pb—1.2%; and Mn—19.0%. 
The organic material was largely fibrous vegetable matter, of which peat was 
the main constituent. 

Subsequent cleaning of the whole tunnel, with high-pressure water jets, 
recovered the major proportion of the deficit (see Fig. 3). The whole of the 
deficit was not recovered because free carbon dioxide, present in the water 
during normal operation, had removed lime from the concrete lining causing 
residual hydrated silica particles to flake off under the action of water jetting. 
The resultant surface was somewhat rougher than the original ‘‘shuttered 
concrete’’ finish. 

Normal generation over the next 18 months produced a still further reduction 
in the head loss due to friction. The two mechanisms responsible for this 
phenomenon are: (1) Grit, normally present in the water, acting as an abrasive 
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CWM RHEIOOL 


Height above sea level -m 


Distance - km 


FIG. 1.—Rheidol Hydroelectric Project: Schematic 


ONE YEARS 
OF POSITION 


FIG. 2.—Biofilm 


to smooth the concrete; and (2) the “‘peat’’ starting to reattach as a slime, 
**filling-in’’ and, thus, smoothing some of the hollows. 

Annual monitoring of the reappearance of the deposit has produced data 
on the rate of the deposition and a correlation between deposit thickness and 
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FIG. 3.—Friction Head Loss as Function of Time 
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FIG. 4.—Friction Head Loss as Function of Biofilm Thickness 


head loss (see Fig. 4). Additionally, we are now able to predict the optimum 
cleaning frequency for maximum generation, some 7 yr-8 yr. 
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WInpD-INDUCED MIXING IN STRATIFIED FLUID 
By C. Kranenburg' 


INTRODUCTION 


In this note a simple model is offered for mixed-layer deepening caused by 
wind action at the free surface of a density-stratified fluid. Solar radiation, 
cooling, advection, and Coriolis effects are not considered. To include an aspect 
of mixing in waters of limited horizontal extent, a streamwise pressure gradient 
is taken into account. 

On certain modeling assumptions, the turbulent energy equation yields a 
dimensionless entrainment rate E, that is inversely proportional to the overall 
Richardson number R,, see the review in Ref. 15. Here E, = u,' aH /dt, R, 
= B/ ua; u, = friction velocity related to the shear stress exerted by the wind; 
H = mixed-layer depth; ¢ = time; and B = total mixed-layer buoyancy defined 


[p.(H) — p,(z)] dz 


in which g = acceleration due to gravity; p = density; p, = reference density; 
and z = vertical downward coordinate (z = 0 at the free surface). The subscript 
o refers to the initial (i.e., before the wind starts to blow) density distribution. 
Of the foregoing modeling assumptions, the central one, namely, that the net 
kinetic-energy input is proportional to ur , has been criticized (14,18). 

A different assumption was put forward by Pollard, et. al (11). These authors 
assumed another overall Richardson number, R,,, to be constant. Here R,,, 
= B/U’ and U = mean mixed-layer velocity; A constant R,,, implies E,«= 
R,'’?. Price (12) re-examined the experiments of Kato and Phillips (6) (linear 
initial density distribution), and Kantha, et. al (5) (two-layer fluid) introducing 
a correction for side-wall friction. Prices’ analysis of these experiments led 
to an almost constant R,,, at larger values of R,, (R,, ~ 0.6, experimental scatter 
in the range 0.5 < R,,, < 0.8), in agreement with Ref. 14. 

An explanation of these findings is given herein starting from a two-layer 
structure of the mixed layer: a homogeneous upper layer in which stratification 
effects are disregarded, and a lower transition layer where near-critical conditions 
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of Tech., Stevinweg 4, 2628 CN Delft, The Netherlands. 
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(i.e., the turbulence is nearly collapsing) prevail. Experimental evidence indicating 
near-critical conditions in this layer (the region including the thermocline) is 
provided by the revealing photographs of Woods (19), by Wyatt (21) and Price 
(13); also see Kundu (8). It is therefore assumed that in the transition layer 


in which u = horizontal velocity component; and overbars indicate variables 
averaged for turbulence. The full-scale measurements reported in Ref. 13 strongly 
support the assumption of a homogeneous upper layer. Apart from the influence 
of the Coriolis force, possible processes contributing to the homogeneity are: 
shear-induced mixing (see Ref. 8 for comparison), wave-induced mixing and 
Langmuir circulations. 


Anatyticat 


The density distribution in the transition layer will be assumed linear, as 
suggested by the measurements reported in Ref. 13. It can be shown using 
elementary calculus of variations that of all possible distributions, the linear 
one yields the smallest entrainment rate (for given depth of the transition layer). 
The relevance of this fact will become clear below. Eq. 3 then shows that 
the velocity profile in the transition layer becomes linear also. 

The balance of mass reads 
l H 
+h) Ap = | [p.(H) — dz 

in which h = depth of homogeneous layer. Adopting the Boussinesq approxima- 
tion, the balance of horizontal momentum becomes 
H ap 


Ox 


dt 


djl 
mans | 


in which ap /ax = streamwise pressure gradient; and u, = additional velocity 
in the homogeneous layer [4 = Au + u, and u,(h,t) = 0]. The undisturbed 
layer is assumed to be deep and nonturbulent. Eqs. 1, 2, 3, and 4 yield 


The positive root chosen in Eq. 6 requires the right side of Eq. 5 to be positive. 
Eliminating Au between Eqs. 5 and 6 gives 


d 
VRP +h | = 
t 


1 
0 


in which R,. = critical value of the gradient Richardson number R, 
ap 
0z 
a 
Au R, H-h 
u, R, H+h 
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in which R = 1/2 (R,/R,.); and » = dummy variable. For simplicity it has 
been assumed in Eq. 7 that u, (z, ¢) can be written as u,(z/h). The pressure-gradient 
parameter y is defined by 


In the case of pure setup without horizontal circulation, y would be unity. 

Although it is possible to include arbitrary time dependence of u, in the 
analysis, the friction velocity will be assumed zero for ¢ < 0 and constant 
for t > 0 from now on. Integrating Eq. 7 then gives [H(0) = 0] 


in which | wn) dn 


0 


Eq. 9 shows that a is related to the entrainment rate in the absence of stratification 
(R =0,h =H) 


=a(l 


Obviously, a may be a function of y. 
Using the turbulence kinetic energy equation, it can be shown that A must 
go to zero for R tending to infinity. Eq. 9 then gives HV R = u, (1 — y) 


t, or upon differentiation (note that in general R depends on H) 


d 
— (HV R) 
dH 


Eqs. 9, 11, and 12 indicate zero entrainment rates for y = 1 (the case y > 
| has already been excluded). The simple shear-flow model does not hold in 
this case, however, since the velocity gradient du/az becomes zero at some 
depth. Diffusion of turbulence kinetic energy should then be considered to obtain 
nonzero entrainment rates. Nevertheless, the very small entrainment rates 
observed by Wu (20), Delft Hydraulics Laboratory (2) and Kit, et al. (7) for 
+y = | and by Imberger, et al. (3) for a reservoir indicate that the model presented 
to some extent correctly reproduces the influence of a streamwise pressure 
gradient. 

Apart from the special cases R = 0 and R-+~, the depth, h, of the homogeneous 
layer has not been specified. A minimal depth, h,,, of the mixed layer can 
be derived from Eq. 9 by taking H as a function of A and letting dH/ah = 
0, while keeping time ¢ fixed. The results are (subscript m refers to the minimum) 


4 
H ap 
ox 
dH 
; u, dt 
u,(l— y)t 
h,, = He (14) 
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Price (13) and Kundu (8) report larger depths A than given by Eq. 14. The 
actual entrainment rates will therefore be larger than minimal. Eq. 13 represents 
an absolute minimum for the layer model elaborated, because of the linear 
density profile assumed for the transition layer. 

In practical situations R is a large number (R, > 50, for instance). Adopting 
a = 0.28 (16), which may be on the large side (1), and R,. = 0.25 (13,17) 
gives a °/R < 0.128. Neglecting a’ in Eq. 13 therefore causes a 6% error 
or less, and differentiating Eq. 13 then again yields Eq. 12. The actual entrainment 
rate may be expected to be less than predicted by Eq. 12, because of the 
restriction at zero R (according to this equation e, would increase undefinitely 
for R tending to zero). Since practically the difference between Eqs. 12 and 
13 is small, Eq. 12 is considered a good approximation. Note that in general 
the entrainment rate is not a simple function of R (or R,). The exceptions to 
this result are considered in the following section. 

Neglecting h with respect to H (as in Eq. 12), the overall Richardson number 
R,, becomes proportional to R,.. For the linear profiles assumed one obtains 


It is noted that R,,, does not depend on the density structure of the undisturbed 
layer and on y(y < 1). Assuming R,. = 0.25 (the theoretical value) to 0.3 (15), 
Eq. 15 gives R,, = 0.5-0.6. These values are close to those indicated in Refs. 
12 and 13. 

Homogeneous Undisturbed Layer.—In this case B and therefore R are constant. 
Eq. 12 then gives 


or E, = 0.71/VR, for R,. = 0.25 and y = 0. On examination, this result 
is found to agree with the experiments of Kantha, et al. (5) up to R, = 200. 
At larger R, the observed entrainment rates drop off. As the authors report, 
this may be caused by molecular effects. Also see Ref. 12 for the influence 
of side-wall friction. 

Initial Density Distribution According to Power Law.—In the case where 


p.(z) = p, + ez® 
in which c, B = constants (8 > —1, cB > 0), Eq. 12 yields 


Comparing Eqs. 16 and 18, the entrainment rate in the case of a linear initial 
density distribution (8 = 1) is only half that for a homogeneous lower layer. 
As noted by Price (12) this result agrees, at least qualitatively, with observations. 
The difference has sometimes been attributed to radiation of internal waves 
in the linearly stratified case (4). The foregoing analysis corroborates that the 
density distributions alone can cause the difference. The linear density profile 
assumed for the transition layer may be unrealistic if 8 = 1. Different profiles 
would yield only slightly different entrainment rates, however. A distribution 
according to a parabola with apex at z = 0, for instance, would give R = 


l-y 
2 
*"B+3VR (1) 
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12/25 (R,/R,.) instead of R = 1/2 (R,/R,.). A comparison between Eq. 18 
(B = 1, y = 0, R,. = 0.25) and the experimental results of Katho and Phillips 
(6) is extrapolated by Price (12) to zero aspect (depth to width) ratio. The 
agreement is satisfactory, although the experimental scatter is considerable. 


Conc.upinc Remarks 


The result of Eqs. 16 and 18, namely that E, is proportional to(1 — y)(R,./R,)'”’, 
can be shown to hold for more general density and velocity profiles by using 
an order of magnitude argument: integrating Eq. 5 gives HU = ua (1 — yt; 
assuming in Eq. 3 that dp/az ~ Ap/H and du/daz ~ U/H gives E, = H/(u,t) 
~ (1 — y)(R,./R,)'’”, and, as was to be expected, R,,, ~ R,.. 

Eq. 15 indicates that R,, < 1, R, being certainly less than 0.5. The flow 
in the mixed layer therefore is supercritical. Entrainment rates in subcritical 
flows will be less than found herein, since subcritical flows seem to be essentially 
more stable. The experiments of Lofquist (9) and most of the experiments 
of Moore and Long (10), for instance, are concerned with subcritical flows. 
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UNCERTAINTY IN STEP-BACKWATER PROFILES” 
Closure by Asok Motayed* and David R. Dawdy,° Members, ASCE 


The writers thank Bradley and Eriksen for setting the record straight by 
pointing out that the step-backwater programs analyzed were ‘‘standard step’’ 
not ‘‘direct step.’’ This indeed eliminated the need for inclusion of an errata 
to this closure. 

To minimize the effects of starting elevation on a study reach, extension 
of study reach downstream, as proposed by the discussers, is a commonly 
used procedure in backwater analysis and was performed by the writers in 
the course of this study. Water surface elevations obtained through such 
techniques were reasonably close to the results of the trial and error method. 
This is expected, since once the effect of the arbitrary starting water surface 
elevation is dissipated over a long reach, the only difference expected to remain 
between the profiles would be that due to inherent differences between the 
computer programs. 

The writers fully agree that the bridge loss computation procedures of different 
programs are likely to be the most important difference between the programs. 
The purpose of this paper, however, was not to enumerate the differences 
between the three programs, per se. It was to examine the inherent uncertainties 
associated with other factors. Therefore, purposefully, the difference in water 
surface elevations caused by different bridge loss computational techniques was 
avoided by selecting a reach without any bridge. The major question raised 
in and the conclusion drawn from the study was that even when all other external 
factors like starting water surface elevations, bridge loss computation, etc. are 
eliminated, step-backwater analyses by the three well-known and commonly 
used computer programs contain uncertainties just because of the different 
computational techniques used in solving the same equation of flow. 


“May, 1979, by Asok Motayed and David R. Dawdy (Proc. Paper 14555). 
*Vice Pres., Sheladia Assoc., Inc., Consulting Engrs., 5711 Sarvis Ave., Riverdale, 
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Open CHANNEL FLow VaryiInG Bep RouGHNeEss* 
Closure by Donald W. Knight” and J. Alasdair Macdonald”® 


The writers would like to thank the various discussers for their valuable 
comments. They agree with Nandana Vittal, M.S. Verma, and K. G. Ranga 
Raju that the division of the cross-sectional area into subareas using isovel 
patterns is somewhat arbitrary. What data there is for three-dimensional flow 
in a corner (15,23,30,31) strongly suggests that the isovel pattern is significantly 
affected by the presence of secondary motions. These secondary flows are 
of course intimately related to the distribution of Reynolds and normal stresses 
in the cross section. Although these secondary flows are relatively small, their 
influence upon the turbulence structure and isovel pattern is large. It therefore 
follows that theoretical models based upon the assumption of one-dimensional 
flow, as most current sidewall correction procedures are, must be less than 
satisfactory in describing a three-dimensional phenomenon. The issue is further 
complicated by the absence of any sound theory accounting for influences of 
geometry (cross-sectional shape) and varying boundary roughness (composite 
roughness). One of the ancillary purposes in publishing the original data was 
to highlight the deficiencies in one of the more popular sidewall correction 
procedures. Fig. 11 adequately shows this point. 

Referring to the correction procedures of Einstein and Vanoni and Brooks, 
the data presented in Fig. 13 do not appear to support the discussers’ statement 
that ‘‘In general, the measured shears are smaller and the difference between 
the measured and computed values does not show any systematic variation 
with B/h.’’ The agreement is certainly poor, but particularly so for low B/h 
values (B/h < 5). There appears to be a systematic deviation in the results 
for low B/h values which arises from the three-dimensionality of the flow, 
and this is not surprising in the light of comments made earlier. Taking a smooth 
rectangular channel as an example, it may be shown that the percentage of 
the shear force carried by the walls, increases markedly at low B/h values. 
(Typically SF, as a percentage varies from 25% at B/h = 5-70% at B/h 
= 1.) It is therefore at these low B/h values that the various sidewall correction 
procedures are best tested. 

The new sidewall correction procedure proposed by Nandana Vittal, M. S. 
Verma, and K. G. Ranga Raju based upon the dimensionless group #B/v appears 
to offer an interesting line of enquiry. The writers are pleased to see that their 
resistance data, when plotted in the form of Fig. 14, appears to fit smoothly 
in with other data. This highlights the value of presenting data in tabular form 
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in published papers. However, as to the universality of this new approach, 
the writers would add a word of caution in the light of earlier comments. It 
is difficult to accept that f,, depends upon “#B/v in the manner shown for all 
possible ratios of bed to wall roughness. Indeed some of the scatter in Fig. 
14 may result from mixing the smooth wall, smooth bed results with the smooth 
wall, rough bed results. Because the legend does not give any details it is 
not possible to draw any firm conclusions. The fact that the Vanoni-Brooks 
method gives comparable results to the new method (Fig. 15) does not in itself 
constitute a rigorous validation unless it is shown that the earlier method is 
itself correct, particularly at low B/h values. However the writers are grateful 
to the discussers for having drawn their attention to this new approach, which 
will be examined further. 

The comments by Peter R. Wormleaton, John Allen, and Panos Hadjipanas 
regarding the Preston tube are helpful and correct. When the flume experiments 
were undertaken, Preston’s original calibration curve was employed, although 
Patel’s criticisms of the original calibration were known. It was decided to 
check the overall shear stress values in the way described in the original paper, 
using both energy slope and velocity readings. The errors shown in Fig. 2(b) 
indicate that the agreement between the various methods was reasonable but 
by no means perfect. The in situ calibration of the Preston Tube was therefore 
accepted as being reasonable in the light of a mean 1,/t, value of 0.976 for 
all 50 experiments, and a variation of 0.927-1.100 by series, and a variation 
of 0.956-1.030 by h/B sets. The fact that the discussers’ results in Fig. 12 
largely agree with those of the writers, indicates that no gross errors were 
introduced by their use of Preston’s original calibration. A more detailed 
examination of boundary shear stress distribution is now being conducted in 
a wind tunnel and Patel’s calibration for the Preston tube has been used in 
this study and will be used in subsequent studies. 

The difficulty in selecting an appropriate bed datum for rippled beds is 
well-known, and details of the writers’ method are given in their earlier paper 
(13). The anomaly referred to in Figs. 4 and 5 arises because of the poor results 
for Experiment 16. Because the hydraulic resistance of the channel did change 
appreciably as the roughness spacing was decreased from one series to the 
next, only every other depth run was carried out for odd numbered series. 
The even numbered series are therefore more indicative of trends. The general 
level of scatter in the shear stress results is not very satisfactory and further 
work is required. A detailed study of momentum transfer in asymmetrical sections 


with varying boundary roughness is currently in progress at the University of 
Birmingham. 


Appenvix.—REFERENCES 


30. Gessner, F. B., and Jones, J. B., “‘On some Aspects of Fully Developed Turbulent 
Flow in Rectangular Channels,”’ Journal of Fluid Mechanics, Vol. 23, Part 4, 1965, 
pp. 689-713. 

31. Melling, A., and Whitelaw, J. H., ‘“‘Turbulent flow in a Rectangular Duct,”’ Journal 
of Fluid Mechanics, Vol. 78, Part 2, 1976, pp. 289-315. 
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Page 1168, Table 1, Col. 1, Experiment 04, Col. 7: Should read 0.198 instead 
of 0.178 

Col. 1, Experiment 23, Col. 3: Should read 9.5 instead of 7.5 
Col. 1, Experiment 32, Col. 13: Should read 131 instead of 181 
Col. 1, Experiment 86, Col. 13: Should read 66.2 instead of 6.2 
Col. 1, Experiment 107, Col. 11: Should read 3.83 instead of 2.83 


SECONDARY FLow AND SHEAR Stress AT River Benps* 
Closure by James C. Bathurst,* Colin R. Thorne,° and Richard D. Hey’ 


The writers thank Gotz for his interesting comments and for the opportunity 
to clarify several points. However, as Gotz does not elaborate on the content 
of the references on which his comments are based, it is uncertain how they 
support his conclusions. Presumably the writers’ measurements and those of 
Gotz show different flow characteristics, possibly because the writers’ measure- 
ments were made in rivers while Gotz’s were apparently made in a laboratory 
flume. The writers’ do not agree that this difference invalidates their conclusions 
regarding: (1) A relationship between some discharge-related term such as 
Reynolds number and shear stress uniformity and the relative strengths of the 
primary and secondary flows; (2) the occurrence of the greatest effect (in relative 
but not absolute terms) of secondary circulation at medium discharge; and (3) 
the dependency of shear stress peak magnitudes and positions on secondary 
currents. These conclusions, which are fully explained in the paper, do not 
seem to be refuted by any evidence in the discusser’s paper and in fact are 
supported by data from a number of recent fluvial studies (6,14,55), most notably 
the full results of the work on the River South Esk in Scotland by J. S. Bridge 
and J. Jarvis (8, 9, and personal communication, 1980). The writers did not 
intend to imply that the secondary flow patterns which they measured are constant 
or steady. Firstly, at a given section the strength of the circulation (particularly 
that of the outer bank cell) fluctuates according to the surging produced by 
the passage of eddies. At the writers’ sites this surging had a period of 5 sec-10 
sec but, as measurements were made for a period of 1 min, the recorded patterns 
are averages for the sections. Secondly, around the bend the pattern varies 
with, among other factors: (1) Cross-sectional shape; (2) outer bank steepness 
{the more likely reason for the difference between Figs. 2(b) and 2(d) noted 


“October, 1979, by James C. Bathurst, Colin R. Thorne, and Richard D. Hey (Proc. 
Paper 14906). 
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by Gotz]; and (3) position along the channel [producing the differences noted 
between Figs. 3(a) and 3(b)]. The considerable changes in pattern around a 
bend and the variations with discharge are described in greater detail elsewhere 
(14, 35, 55, 56 and 57). On the question of how well the results represent 


FIG. 12.—Photograph of Electromagnetic Flowmeter, Wading Rod and Supporting 
Equipment 


bend flow, the writers are uncertain what is meant by this. Natural rivers are 
so full of irregularities and there are so many different flow processes operating 
that, obviously, the precise pattern and scale of the circulation will differ from 
bend to bend. In general, though, the writers’ measurements are consistent 
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with others made at river bends of moderate curvature (e.g., 8, 9, 14, 35, and 
J. S. Bridge and J. Jarvis, personal communication, 1980). 

Concerning field techniques, it is agreed that field measurements are difficult 
to perform accurately but, as secondary velocities in rivers are higher than 
those in laboratory flumes, the measuring instrument does not need to be as 
sensitive. Obviously the electromagnetic flowmeter is too big to be used in 
flumes or close to channel boundaries because of the resulting flow distortion. 
However, it is negligibly small compared with a river cross section and its 
discoid shape (Fig. 12) is designed specifically to minimize distortion effects. 
Further, as it is a rugged instrument, it has a considerable advantage over 
the hot film anemometer. In general, agreement between the patterns measured 
by the electromagnetic flowmeter and the patterns indicated by the Ott C-31 
current meter, movement of tracer material, lines of vortices, and areas of 
upwelling was very good. 

For the measurements the meter was mounted on a specially constructed 
wading rod with a sight and steering arms (Fig. 12) so that it was about 0.3 
m (0.98 ft) in front of the rod. A rope (with tape measure) was set across 
the section, perpendicular to the outer bank, and this was used as the reference 
for longstream and cross-stream directions. [Note that the primary and secondary 
flow directions and velocities were obtained subsequently by a mathematical 
technique (4,56)]. At each vertical along the rope the sight was fixed on the 
same point on the horizon so that the several measurements at the vertical 
had the same direction datum. There was little difficulty in holding the rod 
steady except at high flows when buffeting from eddies caused some movement. 
However, as secondary velocities are higher at high flows, the percentage errors 
in velocity measurement would not have increased in the same proportion as 
the absolute errors. It was not expected that the boat or operator would 
significantly distort the flow around the flowmeter since the meter was held 
well in front of the operator. 

Contrary to Gotz’s implication the writers did not attempt to use the Ott 
C-31 current meter at a height of 0.01 m (0.033 ft) above the bed. The distance 
of 0.01 m was in fact the smallest vertical interval between successive measure- 
ments at a vertical. The point nearest to the bed at which the meter could 
be used was 0.08 m (0.26 ft). Velocities closer to the bed were obtained either 


by an Ott C-1 current meter (as near to the bed as 0.025 m or 0.082 ft) or 
else, as described, by extrapolation. 


55. Bathurst, J. C., “Distribution of Boundary Shear Stress in Rivers,” Adjustments 
of the Fluvial System, D. D. Rhodes and G. P. Williams, eds., Kendall / Hunt Publishing 
Co., Dubuque, Iowa, 1979, pp. 95-116. 

56. Thorne, C. R., Bathurst, J. C., and Hey, R. D., “‘Direktmessungen der Sekundarstré- 
mungen in Flussmaandern’’ (Direct Measurements of Secondary Currents at River 
Meanders), Wasser Wirtschaft, in press. 

57. Thorne, C. R. and Hey, R. D., “‘Direct Measurements of Secondary Currents at 
a River Inflexion Point,’’ Nature, Vol. 280, No. 5719, July, 1979, pp. 226-228. 
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Page 1280, paragraph 3, line 7: Should read ‘‘within 10 mm/s (0.033 fps)’’ 
instead of ‘‘within 10 mm/s (0.33 fps)’’ 

Page 1285, paragraph 2, line 8: Should read “‘i.e., |(O/a,) — @/u,)|/(O/a,),” 
instead of “‘i.e., |(4/u,) — (@/u,)|/(U/ia,),” 

Page 1288, Fig. 8, caption: Should read (Gaged with a Rising Stage) (1 m = 
3.28 ft) instead of [Gaged with a Rising Stage (1 m = 3.28 ft)] 

Page 1289, Table 2, Col. 11 heading: Should read “‘U/(gR)'/*” instead of 
1/29» 

Page 1290, paragraph 2, line 13: Should read ‘‘(Figs. 5, 6, and 7).”’ instead 
of ‘‘(Figs. 1, 5, and 7).”’ 

Page 1294, Ref. 8: Should read ‘‘pp. 303-336.”’ instead of ‘‘pp. 203-336.”’ 

Page 1295, Ref. 37: Should read ‘‘May, 1876,’’ instead of ‘‘May 1976,”’ 

Page 1295, Ref. 38: Should read ‘‘United Kingdom, in 1978,”’ instead of ‘‘United 
Kingdom,”’ 


INCEPTION OF SEDIMENT TRANSPORT” 
Closure by M. Selim Yalin* and Emin Karahan,° Members, ASCE 


The writers would like to thank Posey and Narayanan for their comments 
and complimentary information. 

With regard to the comment of Posey, they would like to point out that 
the purpose of the schematical Fig. 3(b) is merely to show that the flow structure 
in the lower part of the viscous sublayer is the same as that of a laminar 
flow (having the same values of the characteristic parameters at the bed). Fig. 
3(b) is not intended to throw any light on the particle-lift mechanism. On the 
other hand the ‘‘deliberate vagueness’’ of the upper boundary of the viscous 
sublayer in Fig. 3(b) should be sufficient to indicate that the writers do not 
exclude the possibility of the penetration of turbulent fluctuations into the viscous 
sublayer. 

In the writers’ view the schematical analysis presented by Narayanan is 
attractive and informative. The writers have some reasons to believe that with 
the further decrement of X., the common part of the initiation curve (on the 
left of Fig. 5) should tend to become parallel to the X.,-axis. This belief appears 
to be in agreement with the analysis of Narayanan. 

The writers do not think that the transition zone ‘‘is prone to scatter’ due 
to the differences in h/D, as the available data for laminar and hydraulically 
smooth turbulent flows from the same pattern. The divergence between “‘laminar”’ 

“November, 1979, by M. Selim Yalin and Emin Karahan (Proc. Paper 14975). 

°Prof., Dept. of Civ. Engrg., Ellis Hall, Queen’s Univ., Kingston K7L 3N6, Canada. 


°Asst. Prof., Dept. of Civ. Engrg., Technical Univ., Istanbul, Turkey; formerly, 
Post-Doctoral Fellow, Queen’s Univ., Kingston, Ontario, Canada. 
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and ‘‘turbulent’’ initiation curves in Fig. 5 takes place due to the change in 
the regime of a turbulent flow (characterized by X., ~ D ~ k,) and not due 


to the conversion of a laminar flow into a hydraulically smooth) turbulent one 
(characterized by R ~ h). 


Stacep Muttiport Dirrusers* 


Closure by Charles W. Almquist,’ A. M. ASCE, 
and Keith D. Stolzenbach,* M. ASCE 


Brocard has made a useful addition to the theory presented in the original 
paper; data from the Jamestown Nuclear Power Station study presented in his 
Fig. 11 show good agreement with his proposed extension. The writers would 
like to comment briefly on certain aspects of this extension, and to indicate 
an alternative approach which may be applicable when assumptions in the 
proposed extension break down. 

In including the cross-flow entrainment directly into the continuity equation 
for the staged diffuser induced flow, two assumptions have been made; these 
are that the cross-flow entrainment is directly proportional to the cross-flow 
velocity, and that the concept of the staged diffuser as a continuous line source 
of momentum remains valid in the presence of a cross-flow. The former 
assumption has been used with success in the calculation of other jet-type flows 
(23), although the coefficient of proportionality must, in general, be determined 
experimentally and does not have a universal value. For the case of a staged 
diffuser, the value of a. = 1.0 proposed by Brocard appears to be dictated 
by a total flow constraint (see next paragraph) and is not always represented 
of a local dynamic condition as implied in Brocard’s differential entrainment 
relationship (Eq. 29). 

The analytical description of the staged diffuser as a continuous line source 
of momentum is justified in a quiescent ambient by the empirical observation 
that the individual jets of a typical design interact quite strongly, leading in 
fact to an inhibited rate of growth of the diffuser plume along the diffuser 
axis. Even for relatively large jet spacings, the interaction is strong enough 
that the line source assumption still holds. In the presence of a relatively weak 
cross-flow, and so long as the jets are not spaced too far apart, the individual 
jets may be expected to continue to interact quite strongly, and the momentum 
line source concept can be expected to remain approximately correct. In this 


“February, 1980, by Charles W. Almquist and Keith D. Stolzenbach (Proc. Paper 15185). 
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situation, complete entrainment of that portion of the crossflow which passes 
directly over the diffuser would be consistent with the line source nature of 
the diffuser. This is equivalent to a. = 1.0, as Brocard assumes in his proposed 
formulation. 

However, in the presence of a stronger cross-flow, of if the jet spacing is 
relatively large, the individual jets may be deflected enough to prevent interaction, 
each individual jet then behaving much as a free jet in a cross-flow. In this 
case, the concept of the staged diffuser as a line source of momentum can 
no longer be expected to hold, and analysis of the individual jets would be 
more appropriate. This is essentially the approach adopted by Adams and 
Trowbridge (24) and Trowbridge (25) in a detailed investigation of the near-field 
flow of staged diffusers. 

Our conclusion is that, in the presence of a cross-flow, staged diffuser 
performance will fall between two extremes: complete interception of the 
cross-flow, valid for weak cross-flows when the line source concept remains 
valid, and breakdown of the line source concept in the presence of a stronger 
cross-flow, in which case an individual jet analysis is appropriate. Which approach 
is the more applicable must be determined on a case-by-case basis. 

The writers would also like to point out that Adams and Trowbridge (24), 
and Trowbridge (25) have presented comprehensive analyses of the near-field 
of the staged diffuser, including analytical descriptions of the entrainment 
flow-field in both Eulerian and Lagrangian coordinates, and an extension of 
the theory to include sloping bottoms. The Lagrangian analyses form a useful 
basis for evaluating the impact of a staged diffuser discharge on receiving water 


biota; the general approach reported therein may also be applied to other types 
of buoyant discharges. The interested reader is referred to their reports. 


23. Chan, D. T.-L., Lin, J.-T. and Kennedy, J. F., ‘‘Entrainment and Drag Forces of 
Deflected Jets,’’ Journal of the Hydraulics Division, ASCE, Vol. 102, No. HYS, 
Proc. Paper 12141, May, 1976, pp. 615-635. 

24. Adams, E. E., and Trowbridge, J. H., ‘‘Near Field Performance of Staged Diffusers 
in Shallow Water,’’ Energy Laboratory Report No. MIT-EL-79-015, Massachusetts 
Institute of Technology, Cambridge, Mass., Apr., 1979. 

. Trowbridge, J. H., ‘‘Near Field Performance of Staged Diffusers,’’ thesis presented 
to Massachusetts Institute of Technology, at Cambridge, Mass., in 1979, in partial 
fulfillment of the requirements for the degree of Master of Science. 
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Time-DepeNnDENT Stocuastic or Fioops* 


Errata 


The following corrections should be made to the original paper: 


Page 653, Eq. 28, right side: Should read 


n 
0 i=l 


instead of 


| exp - | - 


Page 663, Ref. 5, line 3: Should read ‘‘New York, N.Y., 1972,’’ instead of 
York, N.Y., 1072,” 


Loc Pearson Type 3 Distrisution: 
b 
Metuop or Mixep Moments 


Errata 


The following corrections should be made to the original paper: 


Page 1000, paragraph 4, line 1: Should read ‘‘and c <= y < + o.”’ instead 
of “‘andc = y > + 

Page 1016, Table 8, heading of Col. 7: Should read ‘“‘y,,,.”” instead of “*y,,.,”” 
Page 1017, paragraph 1, line 4: Should read ‘‘By virtue of”’ instead of “‘By 
virture of”’ 


Page 1018, line 1: Should read “‘selection of y, on the basis of Y,”’ instead 
of ‘‘selection y, on the basis of Y,”’ 


“May, 1980, by Michel North (Proc. Paper 15405). 
*June, 1980, by Donthamsetti Veerabhadra Rao (Proc. Paper 15477). 
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Empiricac INVESTIGATION OF Curve NuMBER TECHNIQUE" 
Discussion by Michael Daly,” M. ASCE 


The writer wishes to thank the author for a most informative presentation. 
The author speculates that the reason for the poor fit at the Sonoita Creek 
watershed is that the curve number (CN) technique is inappropriate in the arid 
west. The writer wishes to concur in this speculation as it has been his experience 
that the CN technique provides excessive runoff for large rainfall events. 

An alternate technique developed by the United States Geological Survey 
(USGS) has been employed by the writer which has been found to provide 
very good data (8). This technique uses multiple-regression analysis relating 
flood peaks of 5-yr, 10-yr, 25-yr, and 50-yr recurrence intervals to selected 
physical and climatic basin characteristics. The method is based on data for 
163 sites where flood records have been obtained for 8 yr or more, and on 
the maximum known floods at 439 sites. Using this technique, an estimate 
of the natural-peak flow can be obtained for any desired site within New Mexico 
by using the basin characteristics at the desired site. Copies of the report detailing 
this technique are available from the USGS at the address given in Ref. 8. 
The USGS final report will be available in the summer of 1981. It is being 
prepared by R. P. Thomas. 


8. Scott, A. G., ‘‘Preliminary Flood-Frequency Relations and Summary of Maximum 
Discharges in New Mexico—A Progress Report,’’ Open File Report, U.S. Geological 
Survey, Water Resources Division, Room 115, Federal Building, Santa Fe, N.M. 87501. 


Discussion by Kenneth G. Renard,’ M. ASCE 


The author is to be commended for his effort to present a method for illustrating 
how runoff frequency relationships can be developed from a frequency relation- 
ship for precipitation. 

It would be interesting if the author had commented regarding whether the 
curve numbers he used to develop the information in Figs. 1-5 agreed with 
values suggested in the National Engineering Handbook of the Soil Conservation 
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Service (SCS) (3). Similarly, it would be helpful to note which of the antecedent 
moisture levels was used, as suggested in the handbook. Because the handbook 
suggests that the curve number changes with varying antecedent moisture, it 
seems likely that the curve number associated with the author’s frequency analysis 
would change in response to the exceedence probability. Thus, one might expect 
the curve number for high probability would be lower than it would be for 
low probability. 

That the results of the technique are poor for Sonoita Creek, near Patagonia, 
Ariz., is certainly not surprising to a hydrologist familiar with the hydrologic 
and physiographic characteristics controlling the rainfall-runoff process in the 
region. The spatial distribution of precipitation and the resulting partial area 
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FIG. 6.—Storm near Upper End of Walnut Gulch on July 30, 1966; Each Circle Shows 
Location of Raingage on 57.7-sq mile area 


runoff, along with the reductions in the flow volume (transmission losses) as 
runoff moves from the source area over the normally dry alluvial streambeds, 
dominate the hydrologic response of many semiarid watersheds to individual 
precipitation events. 

To illustrate these phenomena, two runoff events on the Walnut Gulch 
Experimental watershed near Tombstone, Ariz. are used. The watershed is a 
57.7-sq mile (150-km’) tributary of the San Pedro River about 75 mile from 
Sonoita Creek. Mean annual precipitation is about 14 in. (350 mm), which is 
slightly less than that at Sonoita Creek. The hydrologic phenomena, however, 
are very similar in both watersheds. 

Fig. 6 represents a precipitation event concentrated in the upper portion of 
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the watershed. Precipitation was not recorded for this event at the raingage 
in Tombstone. Much of the runoff was lost by transmission losses in the channel 
before reaching the watershed outlet at Flume 1. Such an event shows a typical 
air-mass thunderstorm’s limited areal extent. Runoff measured at the outlet 
of a small 3.18-sq mile (8.24-km”) subwatershed of Walnut Gulch was appreciable. 
By the time this flow had traversed nearly 11.4 mile (18.3 km) of normally 
dry streambed, the flow was significantly less. 

Fig. 7 shows what can happen if the storm is located nearer the watershed 
outlet. For this storm, the opportunity for transmission losses to reduce the 
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FIG. 7.—Storm of August 25, 1968 Was Concentrated in Lower Portion of Watershed; 
Multipeaked Hydrograph Resulted from Runoff in Various Tributaries and from Longer 
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streamflow are much less, and in fact, runoff from different tributaries results 
in a multipeaked hydrograph. If the storm timing was such that all peaks arrived 
at the same time, the flood peak would have been significantly larger. 

In summary, it should be apparent that no single precipitation gage will provide 
an adequate representation of the input to the watershed for a single runoff 
event. Thus, the frequency relationship between individual precipitation events 
and corresponding runoff is nonexistent in ephemeral streams in southwestern 
United States. Given that the raingages at Tombstone were used for these two 
events, it would not have recorded the storm shown in Fig. 6, and would have 
seriously underestimated the precipitation producing the runoff in Fig. 7. 
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CALIBRATION OF Bep-Loap SAMPLERS* 
Discussion by Glendon T. Stevens, Jr.,” M. ASCE 


The authors of this article are to be commended for their display of courage 
in addressing the mysteries of bed-load movements. The same should be said 
for those who, in the past, have probed the recesses of this enigma. They 
have developed an intriguing assortment of samplers and bed load movement 
relationships. Most of these samplers and corresponding bed movement relation- 
ships have been developed with the aid of small-scale physical models. However, 
after having spent extended periods of time at the Waterways Experiment Station 
observing and discussing moveable bed models, it is the writer’s understanding 
that due to distortion in scaling of this type of model, volume of transport 
cannot be predicted. This point is made, however, in light of an understanding 
that this type of model may be utilized to predict movement trends, i.e., location 
of possible scour and deposition sites. 

The primary purpose of this review is two-fold: (1) To inform those concerned 
of the existence of voluminous quantities of prototype data; and (2) to humbly 
probe for answers to some resurgent questions. 

Over the past 10 yr, this writer has been privileged to spend hundreds of 
hours on the Mississippi River collecting data and observing the flow phenomena 
associated with a large alluvial river. During this period this writer collected 
and analyzed a large number of bed samples which were collected with a “‘grab 
type’’ sampler. In addition, with the aid of the latest electronic positioning 
equipment and a recording pathometer this writer collected and observed both 
lateral and longitudinal bed profiles. Longitudinal profiles 1,000 ft in length 
were collected in 100-ft increments across the river at various locations over 
a wide range of flow conditions. 

Analysis of this data reveals that the river bed, at least from the standpoint 
of bed forms, does not display the same activity throughout a lateral section. 
Observation reveals that the height, length, and location of bed activity changes 
under constant flow conditions. In other words, a set of data might reveal 
a segment of the section that sustains sand waves with the remainder of the 
section relatively smooth. One day later the segment that sustained the sand 
waves might be smooth with bed activity in another segment of the section. 
In some cases the bed forms appear laterally yet in others they appear to be 
longitudinal. 

Grab bed samples collected twice daily at the same point in a cross section 
when analyzed show a marked change in D.,. Vertical velocity profile data 
collected simultaneously 5 ft apart when analyzed will produce different coeffi- 


‘October, 1980, by Peter Engel and Y. Lam Lau (Proc. Paper 15725). 
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cients and exponents for a vertical velocity equation. 

And thus, it is that this writer has read published materials dealing with 
sampling, calculating, and predicting sediment transport with highly mixed 
emotions which include interest, chagrin, and wonderment. It appears that each 
new publication sets upon a foundation of fallacies attributable to earlier authors. 
Then, the most recent author tries valiantly to vindicate the historical fallacies 
via implementation of physical model data or mathematical theory. Such procedure 
is highly transparent and equally shocking to the practicing engineer who strives 
to use “‘up-to-date”’ techniques. And each time the practicing engineer is subjected 
to this procedure he has a feeling of having been “‘bitten to death by ducks,” 
since it seems to be a persistent approach utilized by authors today. 

Questions exist on this subject, which, if answered, would clarify an otherwise 
nebulous perspective in regard to the collection and calculation of sediment 
transport. 

First of all, what is the clear-cut definition of bed load? Some authors define 
it as “‘that portion of the sediment load which moves in contact with the river 
bed and possesses a thickness of two-grain diameters.”’ If we are to subscribe 
to this definition, then we have to be concerned about the obvious fact that 
the opening in the basket and sack type samplers appears to be too large! 

Secondly, where does the sampler collect samples? As indicated previously, 
the entire cross section does not appear to be active. There should be places 
in the nonactive zone where bed movement is impending. The presence of 
a sampler in this zone would cause enough turbulence to create motion; thus, 
a sample would be collected which would not be representative of the true 
bed movement. 

Next, how is a sample analyzed? Basket and sack type samplers are fashioned 
like sieves. Therefore, the minimum-sized particle collected is determined by 
the size of the openings in the sampler and the result is that the sample collected 
is not representative of the material transported. 

Another question involves the calculation of bed load; how is this done? 
(Assuming the correct sample has been collected from the proper points within 
the given cross section.) Additionally, what happens when a large sand wave 
(five feet high or greater) passes over a bed load sampler? Also, it is necessary 
to ask whether or not bed load movement has a constant, uniform depth. 

And, finally, two further questions remain unanswered which are worthy 
of consideration: (1) Are total bed load calculations based on the assumption 
that the total depth of bed load has passed through the sampler; and (b) how 
is suspended sediment kept from entering the sampler as it traverses through 
the water? 

Surely these questions are pertinent. If not answerable, then they would tend 
to denigrate the worth of any article addressing the subject and relegate it 
to the category of semantic jousting. 
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SEDIMENTATION IN IRON GATES RESERVOIR ON THE DANUBE* 
Discussion by Nani G. Bhowmik’ 


The author must be congratulated for presenting such a thorough review of 
the sedimentation problems associated with a large reservoir. The size of the 
reservoir, its drainage basin, and other characteristics make it immensely 
important to the well-being of the surrounding countries. 

The author has systematically presented some valuable data as to the suspended 
sediment load, its size distribution, and the deposition of the suspended and 
bed loads in the reservoir. 

The author has made a computation as to the rate of erosion from the catchment 
based on the total sediment load carried by the river. This was shown to be 
3 mm from the entire catchment corresponding to the yearly sediment load 
of 32,200,000 tons in the river. The writer believes that the assumption of 32.2 
x 10° tons of sediment load carried yearly by river is the only sediment that 
is eroded from the catchment is in error. Some of the sediment eroded from 
the catchment does not necessarily reach the main stream. Research conducted 
at the Illinois State Water Survey by Stall and Lee (4) has shown that the 
range of sediment delivery ratios can be anywhere from 20%-50% for the 
watersheds of the State of Illinois. This means that only 20%-50% of the sediment 
eroded from the catchment finally reaches the stream. Therefore, the total weight 
of the eroded sediment from the Danube catchment must be higher than the 
measured sediment load of 32.2 x 10° tons. It would be interesting to find 
out the range of probable sediment delivery ratios for the Danube. 

The author has extrapolated the relationship between river discharge Q and 
sediment load Q, to estimate the sediment load Q, of the Danube of the 133-yr 
period of record. The writer is wondering about the validity of this relationship 
without any constraints. Research conducted at the Illinois State Water Survey 
by Bhowmik, et al. (3) has shown that this type of extrapolation must be done 
with caution. The range of values of sediment load Q, for the same river discharge 
Q can vary by few hundred percent. This is shown for the Wilmington gaging 
station on the Kankakee River in Illinois in Fig. 6. Data shown here are from 
a 12-month period. The confidence intervals of the fitted line are also shown 
in Fig. 6. The catchment area for the station is 1.3 x 10* km’. Similar variabilities 
were also observed for other stations. The writer is very much interested to 
find out about the variabilities observed by the author in his Q@, = Q””/102 
relationship. 

The suspended sediment load analyses shown in Table 5 are also very interesting. 
It appears that most of the suspended sediment load from November-December, 
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1958 at Orsova consisted of silt, clay or fine sand. Data collected from the 
Kankakee River by Bhowmik, et al. (3) indicated that for some stations, during 
flood flows, 60%-80% of the suspended load was sand and the remaining portion 
was silt and clay. However, during low flows the river almost exclusively carried 
silt and clay as suspended load. The writer would like to request the author 
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FIG. 6.—Relationship between Suspended Sediment Load and Water Discharge for 
Kankakee River near Wilmington 


to discuss whether or not such variability was observed for the Danube at 
Orsova or other places. 
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Survey Report of Investigation 98, Urbana, Ill., 1980. 
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